SR_SB_Veg_Analysis

About

Seed rain, seed bank, and vegetational dynamics of north-central Missouri remnant and reconstructed tallgrass prairies

Background:

Author: Katherine Wynne (Wynnekat@msu.edu)

Co-authors: E. Eyerly, L. Sullivan

Created: 14 June 2023

---

Files:

1. Cleaned_Spring_Cover_Data_July_2019.xlsx - contains spring p/a cover data

2. Cleaned_Fall_Cover_Data_September_2019.xlsx - contains fall % cover data

3. Seed_Rain_Data_2019_FINAL.xlsx - contains seed rain data

4. Seed_Bank_2020_Cleaned.xlsx - contains seed bank data

5. SR_Traits_Dataset.xlsx - contains trait data (work in progress)

---

R Version: R 4.2.2

RStudio Version: 2023.03.0+386

Package Version:

Last updated: 16 June 2023

---

—–

Setup

Libraries

### --- Dataset importing and exporting -----

# Import in excel files
library(readxl)

# Export dataframes to excel files

library(writexl)

### ---- Data manipulation and visualization ---- 

# General data cleaning, manipulation, and visualization
library(tidyverse)

## ── Attaching core tidyverse packages ─────────────────── tidyverse 1.3.2.9000 ──
## ✔ dplyr     1.1.0     ✔ readr     2.1.3
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.0     ✔ tibble    3.1.8
## ✔ lubridate 1.9.1     ✔ tidyr     1.3.0
## ✔ purrr     1.0.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors

# 
library(vegan)

## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.6-4

#
library(maditr)

## 
## To drop variable use NULL: let(mtcars, am = NULL) %>% head()
## 
## 
## Attaching package: 'maditr'
## 
## The following objects are masked from 'package:dplyr':
## 
##     between, coalesce, first, last
## 
## The following object is masked from 'package:purrr':
## 
##     transpose
## 
## The following object is masked from 'package:readr':
## 
##     cols

### ---- Model fitting and analysis ---- 

# Mixed-models 

  ## Fits mixed models
    library(lme4)

## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## 
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack

  ## Functions to analyze significance and R^2 components of mixed-models
    library(lmerTest)

## 
## Attaching package: 'lmerTest'
## 
## The following object is masked from 'package:lme4':
## 
##     lmer
## 
## The following object is masked from 'package:stats':
## 
##     step

  ## 
    library(car)

## Loading required package: carData
## 
## Attaching package: 'car'
## 
## The following object is masked from 'package:dplyr':
## 
##     recode
## 
## The following object is masked from 'package:purrr':
## 
##     some

# Summary stats

  ## 
    library(mosaic)

## Registered S3 method overwritten by 'mosaic':
##   method                           from   
##   fortify.SpatialPolygonsDataFrame ggplot2
## 
## The 'mosaic' package masks several functions from core packages in order to add 
## additional features.  The original behavior of these functions should not be affected by this.
## 
## Attaching package: 'mosaic'
## 
## The following objects are masked from 'package:car':
## 
##     deltaMethod, logit
## 
## The following object is masked from 'package:lmerTest':
## 
##     rand
## 
## The following object is masked from 'package:lme4':
## 
##     factorize
## 
## The following object is masked from 'package:Matrix':
## 
##     mean
## 
## The following object is masked from 'package:permute':
## 
##     shuffle
## 
## The following objects are masked from 'package:dplyr':
## 
##     count, do, tally
## 
## The following object is masked from 'package:purrr':
## 
##     cross
## 
## The following object is masked from 'package:ggplot2':
## 
##     stat
## 
## The following objects are masked from 'package:stats':
## 
##     binom.test, cor, cor.test, cov, fivenum, IQR, median, prop.test,
##     quantile, sd, t.test, var
## 
## The following objects are masked from 'package:base':
## 
##     max, mean, min, prod, range, sample, sum

Import Data

setwd("/Users/wynnekat/Desktop/Ph.D. - PFCA and TP Project Datasets")

spring_cover <- read_excel("/Users/wynnekat/Desktop/Ph.D. - PFCA and TP Project Datasets/Vegetative Cover Data 2019/Cleaned_Spring_Cover_Data_July_2019.xlsx")

fall_cover <- read_excel("/Users/wynnekat/Desktop/Ph.D. - PFCA and TP Project Datasets/Vegetative Cover Data 2019/Cleaned_Fall_Cover_Data_September_2019.xlsx")


seed_rain <- read_excel("/Users/wynnekat/Desktop/Ph.D. - PFCA and TP Project Datasets/Seed Rain 2019/Seed_Rain_Data_2019_FINAL.xlsx")

seed_bank <- read_excel("/Users/wynnekat/Desktop/Ph.D. - PFCA and TP Project Datasets/Seed Bank 2020/Seed_Bank_2020_Cleaned.xlsx")

traits <- read_excel("/Users/wynnekat/Desktop/Ph.D. - PFCA and TP Project Datasets/Traits and Species Info/SR_SB_VEG_Traits_Dataset.xlsx")

—–

Methods

Data Collection & Experimental Design

Study Sites

Our study sites were at Tucker Prairie Research Area (38○56’53.6” N, 91○59’40.0” W, Callaway County, MO) and at Prairie Fork Conservation Area (38○58’29.7” N, 91○44’03.3” W, Callaway County, MO). Tucker Prairie is a 59-hectare tract of unplowed North American tallgrass claypan prairie. Less than 0.5 % of intact tallgrass prairie ecosystem (i.e., never-been-plowed) remains in Missouri, and Tucker Prairie represents the last sizable remnant prairie in north central Missouri (Samson & Knopf, 1994). More than 250 species of plants inhabit Tucker Prairie, with representatives from 57 families and over 150 genera (Tropicos, n.d.). Native C4 grasses dominate the landscape including Sorghastrum nutans, Andropogon gerardii, Schizachyrium scoparium, and Sporobolus heterolepis (Drew, 1947; Rabinowitz & Rapp, 1980). Before 1957, Tucker Prairie was used for cattle grazing and occasional haying (Drew, 1947). From 1958 to 2002, Tucker Prairie was burned once every four years in the late winter or early spring (Rabinowitz & Rapp, 1980). Since 2002, Tucker Prairie has been managed on a 5-year burn rotation, where units are burned once in the late winter to early spring (Jan. – Mar.) and again 2-3 years later in the late summer to early fall (Aug. – Oct.).

Prairie Fork Conservation Area (PFCA) is over 450 hectares of former agricultural land being restored to tallgrass prairie and savanna ecosystems (Navarrete-Tindall et al., 2007; Newbold et al., 2019). From 2004 onwards, 16-25 hectares are newly seeded each year with native prairie species collected from Tucker Prairie and other nearby remnant prairies (Newbold et al., 2019). As a result, the plots within PFCA represent a chronosequence of reconstructed tallgrass prairies that are mainly comprised of Tucker Prairie descendants. Prior to seeding, glyphosate-resistant crops (i.e., corn and soybeans) are planted and harvested for at least three years to reduce the existing weed seed bank (Newbold et al., 2019). A seed mix containing an average of 179 species of native forbs, legumes, grasses, rushes, and sedges are broadcast seeded at a rate of 13.4 to 18.2 kg/ha once during the dormant season (Newbold et al., 2019). All mixes are comprised of at least 75% of the same species each year that are representative of nearby remnant prairies such as Tucker Prairie (Newbold et al., 2019). Reconstructions are managed using a 2-4-year burn schedule and annual spraying of invasive species (Newbold et al., 2019). For additional details on PFCA management, see Newbold et al. 2019. To capture changes in seed rain dynamics during the restoration process, we grouped our reconstructed sites into three categories, as defined by Newbold et al. (2019) as being compositionally different. We measured the seed rain in an old reconstruction (seeded in 2004), middle aged reconstruction (seeded in 2012 and 2013), and a young reconstruction (seeded in 2017).

Seed rain sampling

In May 2019, we deployed artificial turf grass seed traps (0.1 x 0.1 m) in Tucker Prairie and at each of the focal PFCA restorations. We used artificial turf grass traps instead of sticky traps due to their durability and resistance to freezing (Molau & Per Mølgaard, 1996; Bass, 2015). Turfgrass traps also do not lose effectiveness over time, like sticky traps, which are prone to contamination by soil, non-seed plant debris, and dead insects, making seed recovery difficult (Arruda et al., 2020). Furthermore, unlike funnel traps, they are easy to install and collect with minimal disturbance to surrounding vegetation and soil and do not fill with water (Schott, 1995). At each site, we randomly established ten, 5 m long transects, placed a minimum of 50 m apart from each other. We then placed traps at 1 m intervals along each transect and affixed them to the ground using ground staples (n = 5 traps per transect, 50 traps per site). We collected and replaced seed traps every 2-weeks during the growing season from May 31st to December 12th, 2019. Hereafter we refer to all captured diaspores as “seeds” despite catching a variety of fruit types (e.g., achene). After collection, all traps from a transect and sampling period were grouped together and sieved through a series of soil sieves (1 mm, 500 mm, 250 mm mesh). From these layers, we then picked out and identified seeds to the lowest taxonomic level possible using a reference collection of species inhabiting our study areas and identification guides (Coons et al., n.d.; A. C. Martin & Barkley, 1961).

Because of their extremely small size (0.3 – 0.5 mm; Yatskievych, 1999), we estimated the number of rush seeds (Juncus sp.) when there were over 200 seeds present in a sample. We first sieved the samples through 180 mm and 150 mm mesh soil sieves. Then we subsampled each sieved layer by calculating the average number of rush seeds per 1 cm2 area (n = 3) and multiplying that average by the total area covered by the sample layer. Lastly, we summed the number of estimated rush seeds per layer to calculate the total number of rush seeds in a sample.

Vegetation sampling

In 2019, at the same transects, we sampled species cover twice, once in the early summer (June - July) and again at peak biomass (September). During the first sampling period, we only recorded species presence. At peak biomass, we also sampled percent aerial species cover in a 1m2 area centered around each seed rain trap. We identified species according to Yatskievych (1999).

Seed bank sampling

In March 2020, before the growing season, we collected soil samples to assess the soil seed bank at each focal prairie. To compare the soil seed bank and the previous year’s seed rain, we collected 5 (10 x 10 x 10 cm3; 1000 cm3) soil cores approximately 0.5 m away from where we captured seed rain at each transect in 2019. We allowed the soil samples to dry before removing all non-seed plant material (e.g., roots and rhizomes) to ensure seedlings only germinated from seeds. We then homogenized all samples collected from a particular transect and subsampled ~1500 cm3 of soil to spread ~1 cm deep over sterile potting soil in plastic germination trays. Starting July 2020, we placed the trays (n = 40) in a greenhouse and watered them when dry. Amongst the 40 trays containing soil samples, we also placed an additional 10 trays containing only sterile potting soil to assess whether any contamination occurred from external sources. Once identifiable, we removed seedlings from trays. After germination ceased in July 2021, we placed all trays into a cold room for ~4 months to replicate conditions necessary to break dormancy for any still dormant seeds. Post-vernalization, we returned the trays to the greenhouse, stirred the soil samples, and monitored for any additional germination. We ended the study one year post-vernalization.

—–

General Dataset Cleaning

# Make the six letter species codes for the spring and fall cover uppercase

spring_cover[[3]] <- toupper(spring_cover[[3]])

fall_cover[[4]] <- toupper(fall_cover[[4]])

# Remove unnecessary columns from datasets

### Spring cover - Only keep Site, Transect, SPP6, Scientific Name

spring_cover <- spring_cover[,-c(5,6)]

### Fall cover - Keep everything

### Seed rain - Only keep Plot, Transect, Sampling Date, Week, SPP6, Number, Unknown

seed_rain <- seed_rain[, -c(7,8)]

### Seed bank - Keep everything


# Change names of columns to better reflect datasets (e.g. seed rain: number -> Number_Seeds)

colnames(seed_bank) <- c("Site", "Transect",  "SPP6", "Number_Seedlings")

colnames(seed_rain) <- c("Site", "Transect",  "Sampling_Date", "Week", "SPP6", "Number_Seeds", "Unknown")


# Change transect PFCA 2 - 10A to just 10 for the cover datasets

spring_cover <- spring_cover %>% 
  mutate_if(is.character, 
                 str_replace_all, pattern = "10A", replacement = "10")

fall_cover <- fall_cover %>%
  mutate_if(is.character, 
                 str_replace_all, pattern = "10a", replacement = "10")


# Make Transect and Plot variables factors for all datasets

## Spring cover

spring_cover$Transect <- as.factor(spring_cover$Transect)

## Fall cover

fall_cover$Transect <- as.factor(fall_cover$Transect)
fall_cover$Plot <- as.factor(fall_cover$Plot)

## Seed Bank

seed_bank$Transect <- as.factor(seed_bank$Transect)

## Seed Rain

seed_rain$Transect <- as.factor(seed_rain$Transect)

Lump everything to transect level data (total number of seeds per transect, total number of seedlings per transect, total amount of fall cover per transect)

####### Clean Spring

# 0.) Make a frequency column 

spring_cover$Cover <- rep(1, nrow(spring_cover)) 

####### Clean Fall

# 0.) Make anything less than 1 = 1

fall_cover <- fall_cover %>% 
  mutate(Cover = ifelse(Cover < 1, 1, Cover))

# 1.) get a sum per species per transect of cover.

fall_sum <- fall_cover %>%
              group_by(Site, Transect, SPP6, Scientific_Name) %>%
              summarize(Cover = sum(Cover))

## `summarise()` has grouped output by 'Site', 'Transect', 'SPP6'. You can
## override using the `.groups` argument.

fall_sum

## # A tibble: 1,215 × 5
## # Groups:   Site, Transect, SPP6 [1,214]
##    Site   Transect SPP6   Scientific_Name          Cover
##    <chr>  <fct>    <chr>  <chr>                    <dbl>
##  1 PFCA 1 1        ACAVIR Acalypha virginica        11.5
##  2 PFCA 1 1        AMBART Ambrosia artemisiifolia    6  
##  3 PFCA 1 1        BARVUL Barbarea vulgaris          3.5
##  4 PFCA 1 1        CHAFAS Chamaecrista fasciculata 388  
##  5 PFCA 1 1        CONCAN Conyza canadensis          1  
##  6 PFCA 1 1        ERISPP Erigeron spp.              2  
##  7 PFCA 1 1        ERIVIL Eriochloa villosa          2  
##  8 PFCA 1 1        EUPSER Eupatorium serotinum       3  
##  9 PFCA 1 1        KUMSTI Kummerowia stipulacea      6  
## 10 PFCA 1 1        LESCAP Lespedeza capitata         1  
## # … with 1,205 more rows

####### Seed Rain

# 1.) get a sum per species per transect of cover.

seed_rain_sum <- seed_rain %>%
              group_by(Site, Transect, SPP6) %>%
              summarize(Number_Seeds = sum(Number_Seeds))

## `summarise()` has grouped output by 'Site', 'Transect'. You can override using
## the `.groups` argument.

# Clean in the same way you did the seed rain analysis

for(i in 1:nrow(seed_rain_sum)){
  if(seed_rain_sum[i,3] == "AGAFAS"){seed_rain_sum[i,3] <- "AGASPP"}
  if(seed_rain_sum[i,3] == "EUPSER"){seed_rain_sum[i,3] <- "EUPSPP"}
  if(seed_rain_sum[i,3] == "GENPUB"){seed_rain_sum[i,3] <- "GENSPP"}
}

####### Seed Bank

# 1.) get a sum per species per transect of cover.

seed_bank_sum <- seed_bank %>%
              group_by(Site, Transect, SPP6) %>%
              summarize(Number_Seedlings = sum(Number_Seedlings))

## `summarise()` has grouped output by 'Site', 'Transect'. You can override using
## the `.groups` argument.

Remove Unknowns from Data sets

####### Remove Spring Unknowns

spring_cover <- filter(spring_cover, !grepl('UNK', SPP6))

# Looks good
list(unique(spring_cover$SPP6))

## [[1]]
##   [1] "ACAVIR" "ACERUB" "ACHMIL" "AGRGIG" "AGRHYE" "ALOCAR" "AMBART" "AMBPSY"
##   [9] "AMOCAN" "ANAMIN" "ANDGER" "ANTNEG" "ASCSYR" "BAPALB" "BAPAUS" "BAPBRA"
##  [17] "BARVUL" "BIDARI" "BLECIL" "BROJAP" "CAMRAD" "CAPBUR" "CARSPP" "CARBIC"
##  [25] "CERSPP" "CHAFAS" "CIRALT" "COMUMB" "CONCAN" "CORLAN" "CORPAL" "CORSPP"
##  [33] "CORTRI" "CROSAG" "CYPACU" "CYPECH" "DALPUR" "DESILL" "DESSES" "DESSPP"
##  [41] "DICLAN" "DIGISC" "ECHPAL" "ELETEN" "ELYVIR" "ERISPP" "ERYYUC" "EUPCOR"
##  [49] "EUPSER" "EUTGYM" "FESARU" "FESPAR" "FRAVIR" "GALAPA" "GALOBT" "GENPUB"
##  [57] "GERCAR" "HELFLE" "HELHEL" "HELMOL" "HELSPP" "HYPPUN" "JUNBRA" "JUNMAR"
##  [65] "JUNSPP" "JUNVIR" "KOEMAC" "KUMSPP" "LACSER" "LEPVIR" "LESCAP" "LESCUN"
##  [73] "LESVIR" "LIAPYC" "LIASQU" "LINSUL" "LOBSPI" "LYSLAN" "LYTALA" "MEDLUP"
##  [81] "MELSPP" "MONFIS" "MUHSPP" "MYOVER" "OENBIE" "OXADIL" "PARINT" "PENDIG"
##  [89] "PLAVIR" "POAPRA" "POLSAN" "POLVER" "POTSIM" "PYCPIL" "PYCTEN" "RANABO"
##  [97] "RATPIN" "RHUCOP" "ROSCAR" "ROSSPP" "RUBSPP" "RUDHIR" "RUDSUB" "RUMCRI"
## [105] "SALAZU" "SCHSCO" "SENMAR" "SETFAB" "SETPAR" "SILANT" "SILINT" "SILLAC"
## [113] "SISSPP" "SOLALT" "SOLCAR" "SOLNEM" "SOLRIG" "SORNUT" "SPHOBT" "SPOHET"
## [121] "STRLEI" "SYMERI" "SYMLAN" "SYMLAT" "SYMNOV" "SYMOBL" "SYMOOL" "SYMORB"
## [129] "SYMPRA" "SYMSPP" "THLARV" "OENSPP" "TRAOHI" "TRIPER" "TRIREP" "TRISPP"
## [137] "ULMSPP" "SCISPP" "VERARV" "VERHEL" "VERSPP" "VIOSAG" "VIOSPP" "VULOCT"
## [145] "ZIZAUR" "SCLTRI"

####### Remove Fall Unknowns

fall_sum <- filter(fall_sum, !grepl('UNK', SPP6))

# Looks good
list(unique(fall_sum$SPP6))

## [[1]]
##   [1] "ACAVIR" "AMBART" "BARVUL" "CHAFAS" "CONCAN" "ERISPP" "ERIVIL" "EUPSER"
##   [9] "KUMSTI" "LESCAP" "PENDIG" "PLAVIR" "RUDHIR" "SETFAB" "SOLALT" "STRLEI"
##  [17] "SYMPIL" "ACHMIL" "CORTRI" "DALPUR" "DIGISC" "EUTGYM" "JUNSPP" "PYCTEN"
##  [25] "RATPIN" "RUDSUB" "SCISPP" "SOLSPE" "SORNUT" "SYMERI" "SYMNOV" "SYMORB"
##  [33] "TRICAM" "BROJAP" "HELMOL" "HYPPUN" "OENBIE" "SOLRIG" "SYMLAN" "AGRHYE"
##  [41] "KUMSTR" "LESCUN" "LIAPYC" "PARINT" "SOLCAR" "TAROFF" "BAPALB" "ERASPE"
##  [49] "JUNVIR" "LESVIR" "SYMPRA" "VERARV" "CERSPP" "ECHCRU" "OXADIL" "PLAMAJ"
##  [57] "SILINT" "SILLAC" "VERSPP" "CORLAN" "CYPECH" "CYPSTR" "EUPCOR" "MEDLUP"
##  [65] "SCHSCO" "BIDARI" "FESPAR" "MONFIS" "RANABO" "SOLNEM" "OENFIL" "ULMPUM"
##  [73] "DESSPP" "ELYCAN" "AGATEN" "BAPBRA" "CARSPP" "CIRALT" "ERYYUC" "GENAND"
##  [81] "KOEMAC" "LACSER" "PYCPIL" "RUMCRI" "SPOHET" "SYMOBL" "VERHEL" "BLECIL"
##  [89] "BOLAST" "CYPACU" "SYMLAT" "SYMOOL" "ECHPAL" "MELSPP" "SYMANO" "SYMLAE"
##  [97] "CORPAL" "ANDGER" "EUPALT" "POAPRA" "ALOCAR" "LIASPP" "LYTALA" "SPOCOM"
## [105] "SPILAC" "PANVIR" "SENMAR" "TRIFLA" "HELSPP" "ELYVIR" "HELHEL" "MYOVER"
## [113] "SALAZU" "TRAOHI" "BAPAUS" "LESPRO" "DESSES" "ROSCAR" "AMOCAN" "LEPDEN"
## [121] "CAMRAD" "RUBSPP" "ACERUB" "ZIZAUR" "DICLAN" "EUPPER" "GALOBT" "LYSLAN"
## [129] "POTSIM" "SETPAR" "SOLMIS" "VIOSAG" "AGAFAS" "ANTNEG" "CROSAG" "LINSUL"
## [137] "RHUGLA" "FRAVIR" "GENPUB" "MUHGLA" "JUNBRA" "ULMSPP" "POLVER" "TRISPP"
## [145] "CORSPP" "POLSAN" "SILANT" "TRIPER" "RHUCOP"

####### Remove Seed Rain Unknowns

# 1.) get a sum per species per transect of cover.


seed_rain_sum <- filter(seed_rain_sum, !grepl('UNK', SPP6))
seed_rain_sum <- filter(seed_rain_sum, !grepl('PANICUM?', SPP6))
seed_rain_sum <- filter(seed_rain_sum, !grepl('NONE', SPP6))

# Looks good
list(unique(seed_rain_sum$SPP6))

## [[1]]
##   [1] "ACAVIR"  "AMBART"  "ANAMIN"  "BOLAST"  "CERSPP"  "CHAFAS"  "CHEALB" 
##   [8] "CONCAN"  "DIGISC"  "ERISPP"  "ERIVIL"  "EUPSPP"  "KUMSPP"  "LEPVIR" 
##  [15] "MEDLUP"  "OXADIL"  "PENDIG"  "PLAVIR"  "PYCTEN"  "RUDHIR"  "SETSPP" 
##  [22] "SOLALT"  "SORNUT"  "SYMSPP"  "THLARV"  "TRIPER"  "VEROSP"  "ALOCAR" 
##  [29] "CORTRI"  "CYPECH"  "DICLAN"  "ERYYUC"  "HYPPUN"  "MONFIS"  "MYOVER" 
##  [36] "OENBIE"  "RATPIN"  "SCIGEO"  "SOLRIG"  "SPHOBT"  "STRLEI"  "TRIFLA" 
##  [43] "ACHMIL"  "BARVUL"  "FESARU"  "GERCAR"  "LESCUN"  "LIAPYC"  "PANCAP" 
##  [50] "SYMNOV"  "BIDARI"  "CARBUS"  "DESSPP"  "ERASPE"  "GALAPA"  "LESVIR" 
##  [57] "POAPRA"  "RUDSUB"  "VERHAS"  "VERSPP"  "ANDGER"  "CARFES"  "ECHCRU" 
##  [64] "EREHIE"  "PLAOCC"  "SILINT"  "AGASPP"  "BROJAP"  "CAPBUR"  "CIRALT" 
##  [71] "CORLAN"  "MELSPP"  "MOLVER"  "SCHSCO"  "SCLTRI"  "CYPSTR"  "FESPAR" 
##  [78] "JUNSPP"  "PYCPIL"  "TAROFF"  "AMATUB"  "LESCAP"  "OENFIL"  "KOEMAC" 
##  [85] "DAUCAR"  "LEUVUL"  "VULOCT"  "CYPACU"  "HORPUS"  "SCIPEN"  "ECHPAL" 
##  [92] "SOLNEM"  "CORPAL"  "LINSUL"  "DIAARM"  "HELMOL"  "PERLON"  "RUMCRI" 
##  [99] "BLECIL"  "LYTALA"  "SPOHET"  "ELESPP1" "SPOCOM"  "CARCEP"  "AGRHYE" 
## [106] "GENSPP"  "HELHEL"  "TRAOHI"  "BAPALB"  "CARANN"  "IPOHED"  "AMOCAN" 
## [113] "SALAZU"  "CARBIC"  "EUPCOR"  "EUTGYM"  "GALOBT"  "RUBSPP"  "SETPAR" 
## [120] "ELESPP"  "LYSLAN"  "POLSPP"  "HYPHIR"  "LYCAME"  "SILANT"  "CARDSP" 
## [127] "CROSAG"  "LOBSPI"  "VIOSAG"

####### Remove Seed Bank Unknowns

seed_bank_sum<- filter(seed_bank_sum, !grepl('UNK', SPP6))

# Looks good
list(sort(unique(seed_bank_sum$SPP6)))

## [[1]]
##  [1] "ABUTHE" "ACAVIR" "ACHMIL" "AGRHYE" "ALOCAR" "AMATUB" "AMBART" "ANDGER"
##  [9] "BARVUL" "CARFES" "CARHIR" "CARPAR" "CERGLO" "CHAFAS" "CIRALT" "CONCAN"
## [17] "CORLAN" "CROSAG" "CYPACU" "CYPSPP" "CYPSTR" "DESSPP" "DICLAN" "DIGISC"
## [25] "DIGSAN" "DIPLAC" "ERASPE" "ERIANN" "ERISTR" "EUPHUM" "EUPSER" "EUTGYM"
## [33] "GERCAR" "HORPUS" "HYPPUN" "IPOHED" "JUNINT" "JUNMAR" "KUMSTI" "KUMSTR"
## [41] "LEPDEN" "LEPVIR" "LESCAP" "LESCUN" "LYSLAN" "MEDLUP" "MELSPP" "MOLVER"
## [49] "OENSPP" "OXADIL" "PANCAP" "PANDIC" "PENDIG" "PERLON" "PLAOCC" "PLAPUS"
## [57] "PLAVIR" "POACHA" "POAPRA" "POPDEL" "POROLE" "PYCPIL" "PYCTEN" "RATPIN"
## [65] "RUDHIR" "RUMCRI" "SCHSCO" "SETFAB" "SETGLA" "SILANT" "SOLALT" "SOLCAR"
## [73] "SORNUT" "SPHOBT" "STRLEI" "SYMLAE" "SYMPIL" "TAROFF" "THLARV" "TOXRAD"
## [81] "TRAOHI" "TRIFLA" "TRIPER" "TRIREP" "VERHAS" "VERPER" "VERTHA"

—–

Merge Spring and Fall cover

# Full join the two cover data sets
cover_merged <- full_join(spring_cover, fall_sum, by = c("Site", "Transect", "SPP6", "Scientific_Name", "Cover")) %>%
                 mutate_all(~replace(., is.na(.), 0))

## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `Transect = (structure(function (..., .x = ..1, .y = ..2, . =
##   ..1) ...`.
## Caused by warning in `[<-.factor`:
## ! invalid factor level, NA generated

# Make sure there is no weirdness in the join

#View(cover_merged) 


### Group by and then select the SPP6 with the max cover between the two datasets

  ### Since Spring only had presence/absence we don't want to inflate the fall cover artificially if we saw that species again

cover_merged_max <- cover_merged   %>%
      group_by(Site, Transect, SPP6, Scientific_Name) %>%
      summarise(Cover=max(Cover))

## `summarise()` has grouped output by 'Site', 'Transect', 'SPP6'. You can
## override using the `.groups` argument.

## Not excatly sure what to do about SYMSPP and KUMSPP (things we couldn't ID in the spring but could in the fall)

    ### Going to drop SYMSPP from the dataset since every plot that had SYMSPP in the spring had identified                    Symphyotrichum in the fall 

cover_merged_max  <- filter(cover_merged_max , !grepl('SYMSPP', SPP6))

    ### Checked KUMSPP only PFCA 2 - 3 and PFCA 2 - 9 did not see Kummerowia again in the fall; will drop KUMSPP for all other plots

          ### Removes all the PFCA 1 KUMSPPs

cover_merged_max <- cover_merged_max[!(cover_merged_max$SPP6%in%"KUMSPP" & cover_merged_max$Site%in% "PFCA 1"),]

         ### Removes all PFCA 2 KUMSPPs that were not in PFCA 2 - 3 and PFCA 2- 9

cover_merged_max <- cover_merged_max[!(cover_merged_max$SPP6%in%"KUMSPP" & cover_merged_max$Site%in% "PFCA 2" & cover_merged_max$Transect%in% c(2,4,5,8)),]


### Check all the species codes

sort(unique(cover_merged_max$SPP6))

##   [1] "ACAVIR" "ACERUB" "ACHMIL" "AGAFAS" "AGATEN" "AGRGIG" "AGRHYE" "ALOCAR"
##   [9] "AMBART" "AMBPSY" "AMOCAN" "ANAMIN" "ANDGER" "ANTNEG" "ASCSYR" "BAPALB"
##  [17] "BAPAUS" "BAPBRA" "BARVUL" "BIDARI" "BLECIL" "BOLAST" "BROJAP" "CAMRAD"
##  [25] "CAPBUR" "CARBIC" "CARSPP" "CERSPP" "CHAFAS" "CIRALT" "COMUMB" "CONCAN"
##  [33] "CORLAN" "CORPAL" "CORSPP" "CORTRI" "CROSAG" "CYPACU" "CYPECH" "CYPSTR"
##  [41] "DALPUR" "DESILL" "DESSES" "DESSPP" "DICLAN" "DIGISC" "ECHCRU" "ECHPAL"
##  [49] "ELETEN" "ELYCAN" "ELYVIR" "ERASPE" "ERISPP" "ERIVIL" "ERYYUC" "EUPALT"
##  [57] "EUPCOR" "EUPPER" "EUPSER" "EUTGYM" "FESARU" "FESPAR" "FRAVIR" "GALAPA"
##  [65] "GALOBT" "GENAND" "GENPUB" "GERCAR" "HELFLE" "HELHEL" "HELMOL" "HELSPP"
##  [73] "HYPPUN" "JUNBRA" "JUNMAR" "JUNSPP" "JUNVIR" "KOEMAC" "KUMSPP" "KUMSTI"
##  [81] "KUMSTR" "LACSER" "LEPDEN" "LEPVIR" "LESCAP" "LESCUN" "LESPRO" "LESVIR"
##  [89] "LIAPYC" "LIASPP" "LIASQU" "LINSUL" "LOBSPI" "LYSLAN" "LYTALA" "MEDLUP"
##  [97] "MELSPP" "MONFIS" "MUHGLA" "MUHSPP" "MYOVER" "OENBIE" "OENFIL" "OENSPP"
## [105] "OXADIL" "PANVIR" "PARINT" "PENDIG" "PLAMAJ" "PLAVIR" "POAPRA" "POLSAN"
## [113] "POLVER" "POTSIM" "PYCPIL" "PYCTEN" "RANABO" "RATPIN" "RHUCOP" "RHUGLA"
## [121] "ROSCAR" "ROSSPP" "RUBSPP" "RUDHIR" "RUDSUB" "RUMCRI" "SALAZU" "SCHSCO"
## [129] "SCISPP" "SCLTRI" "SENMAR" "SETFAB" "SETPAR" "SILANT" "SILINT" "SILLAC"
## [137] "SISSPP" "SOLALT" "SOLCAR" "SOLMIS" "SOLNEM" "SOLRIG" "SOLSPE" "SORNUT"
## [145] "SPHOBT" "SPILAC" "SPOCOM" "SPOHET" "STRLEI" "SYMANO" "SYMERI" "SYMLAE"
## [153] "SYMLAN" "SYMLAT" "SYMNOV" "SYMOBL" "SYMOOL" "SYMORB" "SYMPIL" "SYMPRA"
## [161] "TAROFF" "THLARV" "TRAOHI" "TRICAM" "TRIFLA" "TRIPER" "TRIREP" "TRISPP"
## [169] "ULMPUM" "ULMSPP" "VERARV" "VERHEL" "VERSPP" "VIOSAG" "VIOSPP" "VULOCT"
## [177] "ZIZAUR"

—–

Merge Cover and Seed Rain

### Cover grouping so that it works with the seed rain dataset

# Code to lump certain species together

cover_to_merge_rain <- cover_merged_max

for(i in 1:nrow(cover_to_merge_rain)){
  
  # Group all Symphyotrichum except SYMNOV

  ## SYMPIL
  ## SYMLAN
  ## SYMLAT
  ## SYMERI
  ## SYMANO
  ## SYMLAE
  ## SYMPRA
  ## SYMOBL
  ## SYMOOL

  
  if(cover_to_merge_rain[i,3] == "SYMPIL"){cover_to_merge_rain[i,3] <- "SYMSPP"}
  if(cover_to_merge_rain[i,3] == "SYMLAN"){cover_to_merge_rain[i,3] <- "SYMSPP"}
  if(cover_to_merge_rain[i,3] == "SYMLAT"){cover_to_merge_rain[i,3] <- "SYMSPP"}
  if(cover_to_merge_rain[i,3] == "SYMERI"){cover_to_merge_rain[i,3] <- "SYMSPP"}
  if(cover_to_merge_rain[i,3] == "SYMANO"){cover_to_merge_rain[i,3] <- "SYMSPP"}
  if(cover_to_merge_rain[i,3] == "SYMLAE"){cover_to_merge_rain[i,3] <- "SYMSPP"}
  if(cover_to_merge_rain[i,3] == "SYMPRA"){cover_to_merge_rain[i,3] <- "SYMSPP"}
  if(cover_to_merge_rain[i,3] == "SYMOBL"){cover_to_merge_rain[i,3] <- "SYMSPP"}
  if(cover_to_merge_rain[i,3] == "SYMOOL"){cover_to_merge_rain[i,3] <- "SYMSPP"}
  

# Group all Carex
  ## CARBIC

   if(cover_to_merge_rain[i,3] == "CARBIC"){cover_to_merge_rain[i,3] <- "CARSPP"} 
  
# Group all Kummerowia
  ## KUMSTI
  ## KUMSTR
  
  if(cover_to_merge_rain[i,3] == "KUMSTI"){cover_to_merge_rain[i,3] <- "KUMSPP"} 
  if(cover_to_merge_rain[i,3] == "KUMSTR"){cover_to_merge_rain[i,3] <- "KUMSPP"} 
  
  
# Group all Juncus
  ## JUNMAR
  ## JUNBBRA
  
  if(cover_to_merge_rain[i,3] == "JUNMAR"){cover_to_merge_rain[i,3] <- "JUNSPP"} 
  if(cover_to_merge_rain[i,3] == "JUNBRA"){cover_to_merge_rain[i,3] <- "JUNSPP"} 
  

# Group all Agalinis
  ## AGATEN
  ## AGAFAS
  
  
  if(cover_to_merge_rain[i,3] == "AGATEN"){cover_to_merge_rain[i,3] <- "AGASPP"} 
  if(cover_to_merge_rain[i,3] == "AGAFAS"){cover_to_merge_rain[i,3] <- "AGASPP"} 
  

# Group all Eleocharis
  ## ELETEN 

  if(cover_to_merge_rain[i,3] == "ELETEN"){cover_to_merge_rain[i,3] <- "ELESPP"} 
  
  
# Group all Gentian
  ## GENPUB
  ## GENAND
  
    if(cover_to_merge_rain[i,3] == "GENPUB"){cover_to_merge_rain[i,3] <- "GENSPP"} 
    if(cover_to_merge_rain[i,3] == "GENAND"){cover_to_merge_rain[i,3] <- "GENSPP"} 
  
# Group all Setaria except Setaria parviflora
  ## SETGLA
  ## SETFAB
  
   if(cover_to_merge_rain[i,3] == "SETGLA"){cover_to_merge_rain[i,3] <- "SETSPP"} 
   if(cover_to_merge_rain[i,3] == "SETFAB"){cover_to_merge_rain[i,3] <- "SETSPP"} 
  
  
# Group all Veronica 
  ## VERARV
  
   if(cover_to_merge_rain[i,3] == "VERARV"){cover_to_merge_rain[i,3] <- "VEROSP"} 
  

# Group all Desmodium
  ## DESSES
  
  if(cover_to_merge_rain[i,3] == "DESSES"){cover_to_merge_rain[i,3] <- "DESSPP"} 
  
  
# Group all Eupatorium

  ## EUPALT
  ## EUPSER
  ## EUPPER
  
if(cover_to_merge_rain[i,3] == "EUPALT"){cover_to_merge_rain[i,3] <- "EUPSPP"} 
if(cover_to_merge_rain[i,3] == "EUPSER"){cover_to_merge_rain[i,3] <- "EUPSPP"} 
if(cover_to_merge_rain[i,3] == "EUPPER"){cover_to_merge_rain[i,3] <- "EUPSPP"} 
  
  
  
# Group all Polygala

  ## POLVER
  ## POLSAN
  
if(cover_to_merge_rain[i,3] == "POLVER"){cover_to_merge_rain[i,3] <- "POLSPP"} 
if(cover_to_merge_rain[i,3] == "POLSAN"){cover_to_merge_rain[i,3] <- "POLSPP"} 

  
}  

### Check that it changed the codes

sort(unique(cover_to_merge_rain$SPP6))

##   [1] "ACAVIR" "ACERUB" "ACHMIL" "AGASPP" "AGRGIG" "AGRHYE" "ALOCAR" "AMBART"
##   [9] "AMBPSY" "AMOCAN" "ANAMIN" "ANDGER" "ANTNEG" "ASCSYR" "BAPALB" "BAPAUS"
##  [17] "BAPBRA" "BARVUL" "BIDARI" "BLECIL" "BOLAST" "BROJAP" "CAMRAD" "CAPBUR"
##  [25] "CARSPP" "CERSPP" "CHAFAS" "CIRALT" "COMUMB" "CONCAN" "CORLAN" "CORPAL"
##  [33] "CORSPP" "CORTRI" "CROSAG" "CYPACU" "CYPECH" "CYPSTR" "DALPUR" "DESILL"
##  [41] "DESSPP" "DICLAN" "DIGISC" "ECHCRU" "ECHPAL" "ELESPP" "ELYCAN" "ELYVIR"
##  [49] "ERASPE" "ERISPP" "ERIVIL" "ERYYUC" "EUPCOR" "EUPSPP" "EUTGYM" "FESARU"
##  [57] "FESPAR" "FRAVIR" "GALAPA" "GALOBT" "GENSPP" "GERCAR" "HELFLE" "HELHEL"
##  [65] "HELMOL" "HELSPP" "HYPPUN" "JUNSPP" "JUNVIR" "KOEMAC" "KUMSPP" "LACSER"
##  [73] "LEPDEN" "LEPVIR" "LESCAP" "LESCUN" "LESPRO" "LESVIR" "LIAPYC" "LIASPP"
##  [81] "LIASQU" "LINSUL" "LOBSPI" "LYSLAN" "LYTALA" "MEDLUP" "MELSPP" "MONFIS"
##  [89] "MUHGLA" "MUHSPP" "MYOVER" "OENBIE" "OENFIL" "OENSPP" "OXADIL" "PANVIR"
##  [97] "PARINT" "PENDIG" "PLAMAJ" "PLAVIR" "POAPRA" "POLSPP" "POTSIM" "PYCPIL"
## [105] "PYCTEN" "RANABO" "RATPIN" "RHUCOP" "RHUGLA" "ROSCAR" "ROSSPP" "RUBSPP"
## [113] "RUDHIR" "RUDSUB" "RUMCRI" "SALAZU" "SCHSCO" "SCISPP" "SCLTRI" "SENMAR"
## [121] "SETPAR" "SETSPP" "SILANT" "SILINT" "SILLAC" "SISSPP" "SOLALT" "SOLCAR"
## [129] "SOLMIS" "SOLNEM" "SOLRIG" "SOLSPE" "SORNUT" "SPHOBT" "SPILAC" "SPOCOM"
## [137] "SPOHET" "STRLEI" "SYMNOV" "SYMORB" "SYMSPP" "TAROFF" "THLARV" "TRAOHI"
## [145] "TRICAM" "TRIFLA" "TRIPER" "TRIREP" "TRISPP" "ULMPUM" "ULMSPP" "VERHEL"
## [153] "VEROSP" "VERSPP" "VIOSAG" "VIOSPP" "VULOCT" "ZIZAUR"

### Sum everything by species by transect again 


cover_to_merge_rain <- cover_to_merge_rain %>%
              group_by(Site, Transect, SPP6, Scientific_Name) %>%
              summarize(Cover = sum(Cover))

## `summarise()` has grouped output by 'Site', 'Transect', 'SPP6'. You can
## override using the `.groups` argument.

### Seed grouping so that it works with the cover dataset

seed_rain_names <- left_join(seed_rain_sum, traits)

## Joining with `by = join_by(SPP6)`

rain_to_merge_cover <- seed_rain_names[,c(1,2,3,4,5)]

for(i in 1:nrow(rain_to_merge_cover)){
  
# Group all Carex
  ## CARBUS
  ## CARBIC
  ## CARANN
  ## CARFES
  ## CARCEP
  
  
  if(rain_to_merge_cover[i,3] == "CARBUS"){rain_to_merge_cover[i,3] <- "CARSPP"}
  if(rain_to_merge_cover[i,3] == "CARBIC"){rain_to_merge_cover[i,3] <- "CARSPP"}
  if(rain_to_merge_cover[i,3] == "CARANN"){rain_to_merge_cover[i,3] <- "CARSPP"}
  if(rain_to_merge_cover[i,3] == "CARFES"){rain_to_merge_cover[i,3] <- "CARSPP"}
  if(rain_to_merge_cover[i,3] == "CARCEP"){rain_to_merge_cover[i,3] <- "CARSPP"}
  
 # Group all Eleocharis
  ## ELESPP1
 
  if(rain_to_merge_cover[i,3] == "ELESPP1"){rain_to_merge_cover[i,3] <- "ELESPP"}
  
  
   # Group all Scirpus
  ## ELESPP1
 
  if(rain_to_merge_cover[i,3] == "SCIGEO"){rain_to_merge_cover[i,3] <- "SCISPP"}
  if(rain_to_merge_cover[i,3] == "SCIPEN"){rain_to_merge_cover[i,3] <- "SCISPP"}


}


### Sum everything by species by transect again 


rain_to_merge_cover <- rain_to_merge_cover %>%
              group_by(Site, Transect, SPP6) %>%
              summarize(Number_Seeds = sum(Number_Seeds))

## `summarise()` has grouped output by 'Site', 'Transect'. You can override using
## the `.groups` argument.

# Looks good
sort(unique(rain_to_merge_cover$SPP6))

##   [1] "ACAVIR" "ACHMIL" "AGASPP" "AGRHYE" "ALOCAR" "AMATUB" "AMBART" "AMOCAN"
##   [9] "ANAMIN" "ANDGER" "BAPALB" "BARVUL" "BIDARI" "BLECIL" "BOLAST" "BROJAP"
##  [17] "CAPBUR" "CARDSP" "CARSPP" "CERSPP" "CHAFAS" "CHEALB" "CIRALT" "CONCAN"
##  [25] "CORLAN" "CORPAL" "CORTRI" "CROSAG" "CYPACU" "CYPECH" "CYPSTR" "DAUCAR"
##  [33] "DESSPP" "DIAARM" "DICLAN" "DIGISC" "ECHCRU" "ECHPAL" "ELESPP" "ERASPE"
##  [41] "EREHIE" "ERISPP" "ERIVIL" "ERYYUC" "EUPCOR" "EUPSPP" "EUTGYM" "FESARU"
##  [49] "FESPAR" "GALAPA" "GALOBT" "GENSPP" "GERCAR" "HELHEL" "HELMOL" "HORPUS"
##  [57] "HYPHIR" "HYPPUN" "IPOHED" "JUNSPP" "KOEMAC" "KUMSPP" "LEPVIR" "LESCAP"
##  [65] "LESCUN" "LESVIR" "LEUVUL" "LIAPYC" "LINSUL" "LOBSPI" "LYCAME" "LYSLAN"
##  [73] "LYTALA" "MEDLUP" "MELSPP" "MOLVER" "MONFIS" "MYOVER" "OENBIE" "OENFIL"
##  [81] "OXADIL" "PANCAP" "PENDIG" "PERLON" "PLAOCC" "PLAVIR" "POAPRA" "POLSPP"
##  [89] "PYCPIL" "PYCTEN" "RATPIN" "RUBSPP" "RUDHIR" "RUDSUB" "RUMCRI" "SALAZU"
##  [97] "SCHSCO" "SCISPP" "SCLTRI" "SETPAR" "SETSPP" "SILANT" "SILINT" "SOLALT"
## [105] "SOLNEM" "SOLRIG" "SORNUT" "SPHOBT" "SPOCOM" "SPOHET" "STRLEI" "SYMNOV"
## [113] "SYMSPP" "TAROFF" "THLARV" "TRAOHI" "TRIFLA" "TRIPER" "VERHAS" "VEROSP"
## [121] "VERSPP" "VIOSAG" "VULOCT"

# Drop the scientific names they are causing trouble with the join

cover_to_merge_rain <- cover_to_merge_rain[, -4]


cover_rain <- full_join(cover_to_merge_rain, rain_to_merge_cover) %>%
                 mutate_all(~replace(., is.na(.), 0))

## Joining with `by = join_by(Site, Transect, SPP6)`
## `mutate_all()` ignored the following grouping variables:

# Drop TP transect 7 because of missing data

cover_rain <- cover_rain[!(cover_rain$Site%in%"TP" & cover_rain$Transect%in% "7"),]

—–

Merge Cover and Seed Bank

### Cover grouping so that it works with the seed bank dataset

# Code to lump certain species together

cover_to_merge_bank <- cover_merged_max

for(i in 1:nrow(cover_to_merge_bank)){
  
  # Group all Desmodium together

  ## DESSES

  if(cover_to_merge_bank[i,3] == "DESSES"){cover_to_merge_bank[i,3] <- "DESSPP"}

  # Group all Cyperus together

  ## CYPSTR
  ## CYPECH

  if(cover_to_merge_bank[i,3] == "CYPSTR"){cover_to_merge_bank[i,3] <- "CYPSPP"}
  if(cover_to_merge_bank[i,3] == "CYPECH"){cover_to_merge_bank[i,3] <- "CYPSPP"}

  # Group all Oenothera together

  ## OENFIL
  ## OENBIE
  
  if(cover_to_merge_bank[i,3] == "OENFIL"){cover_to_merge_bank[i,3] <- "OENSPP"}
  if(cover_to_merge_bank[i,3] == "OENBIE"){cover_to_merge_bank[i,3] <- "OENSPP"}



}

cover_to_merge_bank <- cover_to_merge_bank %>% 
  group_by(Site, Transect, SPP6) %>% 
  summarize(Cover = sum(Cover))

## `summarise()` has grouped output by 'Site', 'Transect'. You can override using
## the `.groups` argument.

### Seed Bank grouping so that it works with the cover dataset

bank_to_merge_cover <- seed_bank_sum

# Code to lump certain species together

for(i in 1:nrow(bank_to_merge_cover)){
  
  # Group all Cerastium together

  ## CERGLO

  if(bank_to_merge_cover[i,3] == "CERGLO"){bank_to_merge_cover[i,3] <- "CERSPP"}

  # Group all Cyperus together

  ## DESSES

  if(bank_to_merge_cover[i,3] == "CYPSTR"){bank_to_merge_cover[i,3] <- "CYPSPP"}
  
  # Group all the Carex together

  if(bank_to_merge_cover[i,3] == "CARFES"){bank_to_merge_cover[i,3] <- "CARSPP"}
  
  # Group all Erigeron together

  ## ERIANN
  ## ERISTR

  if(bank_to_merge_cover[i,3] == "ERIANN"){bank_to_merge_cover[i,3] <- "ERISPP"}
  if(bank_to_merge_cover[i,3] == "ERISTR"){bank_to_merge_cover[i,3] <- "ERISPP"}

  # Change Juncus interior to Juncus sp. (since we were able to ID JUNMAR in the cover but not JUNINT)
  
  if(bank_to_merge_cover[i,3] == "JUNINT"){bank_to_merge_cover[i,3] <- "JUNSPP"}
  
}

bank_to_merge_cover <- bank_to_merge_cover %>% 
  group_by(Site, Transect, SPP6) %>% 
  summarize(Number_Seedlings = sum(Number_Seedlings))

## `summarise()` has grouped output by 'Site', 'Transect'. You can override using
## the `.groups` argument.

# Drop the scientific names they are causing trouble with the join



cover_bank <- full_join(cover_to_merge_bank, bank_to_merge_cover) %>%
                 mutate_all(~replace(., is.na(.), 0))

## Joining with `by = join_by(Site, Transect, SPP6)`
## `mutate_all()` ignored the following grouping variables:

# Drop TP transect 7 because of missing fall cover data

cover_bank <- cover_bank[!(cover_bank$Site%in%"TP" & cover_bank$Transect%in% "7"),]

—–

Merge Seed Rain and Seed Bank

### Seed Bank grouping so that it works with the cover dataset

bank_to_merge_rain <- seed_bank_sum

# Code to lump certain species together

for(i in 1:nrow(bank_to_merge_rain)){
  
  # Group all Cerastium together

  ## CERGLO

  if(bank_to_merge_rain[i,3] == "CERGLO"){bank_to_merge_rain[i,3] <- "CERSPP"}

  # Group all Cyperus together

  ## CYPSTR
  
  if(bank_to_merge_rain[i,3] == "CYPSTR"){bank_to_merge_rain[i,3] <- "CYPSPP"}
  
  # Group all Erigeron together

  ## ERIANN
  ## ERISTR

  if(bank_to_merge_rain[i,3] == "ERIANN"){bank_to_merge_rain[i,3] <- "ERISPP"}
  if(bank_to_merge_rain[i,3] == "ERISTR"){bank_to_merge_rain[i,3] <- "ERISPP"}

  # Group all the Juncus together
  
  ## JUNINT
  ## JUNMAR
  
  if(bank_to_merge_rain[i,3] == "JUNINT"){bank_to_merge_rain[i,3] <- "JUNSPP"}
  if(bank_to_merge_rain[i,3] == "JUNMAR"){bank_to_merge_rain[i,3] <- "JUNSPP"}
  
  # Group all the Setaria together
  
  if(bank_to_merge_rain[i,3] == "SETFAB"){bank_to_merge_rain[i,3] <- "SETSPP"}
  if(bank_to_merge_rain[i,3] == "SETGLA"){bank_to_merge_rain[i,3] <- "SETSPP"}
  
  # Group all the Symphyotrichum together
  
  ## SYMPIL
  ## SYMLAE
  
  if(bank_to_merge_rain[i,3] == "SYMPIL"){bank_to_merge_rain[i,3] <- "SYMSPP"}
  if(bank_to_merge_rain[i,3] == "SYMLAE"){bank_to_merge_rain[i,3] <- "SYMSPP"}
  if(bank_to_merge_rain[i,3] == "SIMPIL"){bank_to_merge_rain[i,3] <- "SYMSPP"}
  
   # Group all the Kummerowia together
  
  ## KUMSTR
  ## KUMSTI
  
  if(bank_to_merge_rain[i,3] == "KUMSTR"){bank_to_merge_rain[i,3] <- "KUMSPP"}
  if(bank_to_merge_rain[i,3] == "KUMSTI"){bank_to_merge_rain[i,3] <- "KUMSPP"}
  
   # Group all the Veronica together
  
  ## VERPER

  if(bank_to_merge_rain[i,3] == "VERPER"){bank_to_merge_rain[i,3] <- "VEROSP"}

  # Group all the Eupatorium together
  
  ## EUPSER

  if(bank_to_merge_rain[i,3] == "EUPSER"){bank_to_merge_rain[i,3] <- "EUPSPP"}
  
  # Merge all the Cardamine together
  
  ## CARHIR
  ## CARPAR

  if(bank_to_merge_rain[i,3] == "CARHIR"){bank_to_merge_rain[i,3] <- "CARDSP"}
  if(bank_to_merge_rain[i,3] == "CARPAR"){bank_to_merge_rain[i,3] <- "CARDSP"}
}

bank_to_merge_rain <- bank_to_merge_rain %>% 
  group_by(Site, Transect, SPP6) %>% 
  summarize(Number_Seedlings = sum(Number_Seedlings))

## `summarise()` has grouped output by 'Site', 'Transect'. You can override using
## the `.groups` argument.

rain_to_merge_bank <- seed_rain_sum

for(i in 1:nrow(rain_to_merge_bank)){

  # Group all Cyperus together

  ## CYPSTR
  ## CYPECH

  if(rain_to_merge_bank[i,3] == "CYPSTR"){rain_to_merge_bank[i,3] <- "CYPSPP"}
  if(rain_to_merge_bank[i,3] == "CYPECH"){rain_to_merge_bank[i,3] <- "CYPSPP"}

  # Group all Oenothera together

  ## OENFIL
  ## OENBIE
  
  if(rain_to_merge_bank[i,3] == "OENFIL"){rain_to_merge_bank[i,3] <- "OENSPP"}
  if(rain_to_merge_bank[i,3] == "OENBIE"){rain_to_merge_bank[i,3] <- "OENSPP"}

}


rain_to_merge_bank <- rain_to_merge_bank %>% 
  group_by(Site, Transect, SPP6) %>% 
  summarize(Number_Seeds = sum(Number_Seeds))

## `summarise()` has grouped output by 'Site', 'Transect'. You can override using
## the `.groups` argument.

rain_bank <- full_join(bank_to_merge_rain, rain_to_merge_bank) %>%
                 mutate_all(~replace(., is.na(.), 0))

## Joining with `by = join_by(Site, Transect, SPP6)`
## `mutate_all()` ignored the following grouping variables:

—-

Seed Rain and Cover Exploratory Analysis

cover_rain_sum <- cover_rain %>% 
  group_by(SPP6) %>% 
  summarize(tot.cover = sum(Cover), tot.seeds = sum(Number_Seeds))

Wasn’t in the cover at all but we caught seeds

• Amaranthus tuberculatus
• Cardamine spp.
• Chenopodium album
• Daucus carota
• Dianthus armeria
• Erechtites hieracifolia
• Hordeum pusillum
• Ipomoea hederacea
• Leucanthemum vulgare
• Mollugo verticillata
• Panicum capillare
• Persicaria longiseta
• Platanus occidentalis
• Verbena hastata

Didn’t catch seeds but was in the cover

Notable (> 10% cover observed)

• Baptisia australis
• Potentilla simplex
• Silphium laciniatum
• Rosa carolina
• Fragraria virginica
• Campsis radicans
• Parthenium integrifolium
• Rhus copallinum
• Solanum carolinense
• Helianthus spp (the pressed one)
• Baptisia bracteata
• Antennaria neglecta

Minor ( 1% < cover < 10% observed)

• Elymus virginica
• Rhus glabra
• Trifolium campestre
• Verbesina helianthoides
• Solidago speciosa
• Lespedeza procombens
• Liatris spp.
• Ranunculus abortivus
• Juniperus virginiana
• Trifolium spp.
• Ulmus spp.
• Agrostis gigantea
• Lactuca serriola
• Symphoricarpos orbiculatus
• Ulmus pumila
• Acer rubrum
• Comandra umbellata
• Dalea purpurea
• Muhlenbergia spp.
• Plantago major
• Rosa spp.
• Spiranthes lacera
• Trifolium repens
• Viola spp.
• Zizea aurea

Trace (cover 1% observed)

• Ambrosia psilostachya
• Asclepias syriaca
• Desmanthus illinoensis
• Elymus canadense
• Helenium flexuosum
• Lepidium densiflorum
• Liatris squarrosa
• Muhlenbergia glabrifolia
• Oenothera spp. (not biennis)
• Panicum virgatum
• Senna marilandica
• Sisyrinchium spp.
• Solidago missouriensis.

Dissimilarity analysis

## # A tibble: 2,750 × 5
## # Groups:   Site, Transect, SPP6 [1,809]
##    Site   Transect SPP6   Abundance Type       
##    <chr>  <fct>    <chr>      <dbl> <chr>      
##  1 PFCA 1 1        ACAVIR      11.5 Aboveground
##  2 PFCA 1 1        ALOCAR       1   Aboveground
##  3 PFCA 1 1        AMBART       6   Aboveground
##  4 PFCA 1 1        BARVUL       3.5 Aboveground
##  5 PFCA 1 1        BROJAP       1   Aboveground
##  6 PFCA 1 1        CERSPP       1   Aboveground
##  7 PFCA 1 1        CHAFAS     388   Aboveground
##  8 PFCA 1 1        CONCAN       1   Aboveground
##  9 PFCA 1 1        DESILL       1   Aboveground
## 10 PFCA 1 1        ERISPP       2   Aboveground
## # … with 2,740 more rows

# Code obtained and modified from Lauren and Anu's seed rain vs. vegetation paper


##### run dissimilarity loops
Cover_rain_dis_results <-matrix(nrow=0,ncol=7)
sites <-as.vector(unique(cover_rain_dis$Site))


for(s in 1:length(sites)){
  temp <-subset(cover_rain_dis, Site==sites[s])
  temp$SPP6 <- factor(as.character(temp$SPP6))
  
  transects <-as.vector(unique(temp$Transect))
  
  #for(b in 1:length(blocks)){
    #temp_b<-subset(temp, block==blocks[b])

    #plots<-as.vector(unique(temp_b$plot))
   
      for(t in 1:length(transects)){ 
        
        temp_t <- subset(temp, Transect==transects[t])

        #wide form - with a row for above and belowground
        plot_cast<-dcast(temp_t[,-c(1:2)],Type ~ SPP6, value.var="Abundance", fun.aggregate=mean,fill=0)
        
        #relative abundances
        trans_tot <- cbind(plot_cast$Type,decostand(plot_cast[,-1],"tot"))
        names(trans_tot)[1]<-"trt"
        
        trans_pa <- cbind(plot_cast$Type,decostand(plot_cast[,-1],"pa"))
        names(trans_pa)[1]<-"trt"
        
        trans_int <- cbind(plot_cast$Type,trunc(plot_cast[,-1])) #integers only
        names(trans_int)[1]<-"trt"
        
        trans_int_norm <- cbind(plot_cast$Type,trunc(decostand(plot_cast[,-1], "normalize")*100)) #integers only and normalize
        names(trans_int)[1]<-"trt"
        
        distance_bray<-vegdist(trans_tot[,-1], method = "bray")
        distance_jaccard <- vegdist(trans_pa[,-1], method = "jaccard")
        distance_kulczynski<-vegdist(trans_tot[,-1], method = "kulczynski")
        distance_cao<-vegdist(trans_int_norm[,-1], method = "cao")
        distance_morisita<-vegdist(trans_int[,-1], method = "morisita")
        

        new.row<-c(sites[s], transects[t],distance_bray[1],distance_jaccard[1],
                   distance_cao[1], distance_morisita[1], distance_kulczynski[1])
        Cover_rain_dis_results <-rbind(Cover_rain_dis_results, new.row)
    }
  print(s/length(sites))
}

## [1] 0.25
## [1] 0.5
## [1] 0.75
## [1] 1

row.names(Cover_rain_dis_results)<-NULL
Cover_rain_dis_results<- data.frame(Cover_rain_dis_results)

names(Cover_rain_dis_results) <- c("Site", "Transect", "dissimilarity_bray", "dissimilarity_jaccard",
                    "dissimilarity_cao", "dissimilarity_morisita", "dissimilarity_kulczynski")

## 
## Call:
## lm(formula = dissimilarity_bray ~ Site, data = Cover_rain_dis_results)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.267522 -0.063391 -0.002319  0.088438  0.190470 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.46881    0.03621  12.949 6.58e-15 ***
## SitePFCA 2   0.16217    0.05120   3.167 0.003186 ** 
## SitePFCA 3   0.21871    0.05120   4.272 0.000141 ***
## SiteTP       0.21184    0.05261   4.027 0.000289 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1145 on 35 degrees of freedom
## Multiple R-squared:  0.4014, Adjusted R-squared:  0.3501 
## F-statistic: 7.823 on 3 and 35 DF,  p-value: 0.0003994

## Anova Table (Type III tests)
## 
## Response: dissimilarity_bray
##              Sum Sq Df  F value    Pr(>F)    
## (Intercept) 2.19783  1 167.6676 6.577e-15 ***
## Site        0.30764  3   7.8231 0.0003994 ***
## Residuals   0.45879 35                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Fig. - Bray-Curtis dissimilarity

## 95% confidence intervals 

emmeans::emmeans(Cover_rain.mod.bray, "Site", type = "response")

##  Site   emmean     SE df lower.CL upper.CL
##  PFCA 1  0.469 0.0362 35    0.395    0.542
##  PFCA 2  0.631 0.0362 35    0.557    0.704
##  PFCA 3  0.688 0.0362 35    0.614    0.761
##  TP      0.681 0.0382 35    0.603    0.758
## 
## Confidence level used: 0.95

## Make a new dataframe to hold the model estimates and 95% confidence interval bounds

response <- c(0.472, 0.628, 0.687, 0.681)
lower <- c(0.396, 0.552,  0.611, 0.601) 
upper <- c(0.548, 0.704, 0.763,  0.761)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")


Cover_rain_dissimilarity <- data.frame(site, response, lower, upper)

Fig. - Jaccard dissimilarity figure

## 95% confidence intervals 

emmeans::emmeans(Cover_rain.mod.jac, "Site", type = "response")

##  Site   emmean     SE df lower.CL upper.CL
##  PFCA 1  0.480 0.0193 35    0.440    0.519
##  PFCA 2  0.462 0.0193 35    0.423    0.502
##  PFCA 3  0.589 0.0193 35    0.550    0.628
##  TP      0.542 0.0204 35    0.500    0.583
## 
## Confidence level used: 0.95

## Make a new dataframe to hold the model estimates and 95% confidence interval bounds

response <- c(0.480, 0.462, 0.589, 0.542)
lower <- c(0.440, .423,  0.550, 0.500) 
upper <- c(0.519, 0.502, 0.628,  0.583)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")


Cover_rain_jac_dissimilarity <- data.frame(site, response, lower, upper)

Immigrants

cover_rain_dis_wide <- cover_rain_dis %>% 
  spread(SPP6, Abundance) %>% 
  mutate_all(~replace(., is.na(.), 0)) 

## `mutate_all()` ignored the following grouping variables:
## • Columns `Site`, `Transect`
## ℹ Use `mutate_at(df, vars(-group_cols()), myoperation)` to silence the message.

cover_rain_dis_long <- cover_rain_dis_wide %>% 
  gather(key = "SPP6", value = "Abundance", 4:173) %>% 
  spread(key = Type, Abundance) 

# Get rid of columns that both have zeros
cover_rain_dis_long_unique <- cover_rain_dis_long[!(cover_rain_dis_long$Aboveground == 0 & cover_rain_dis_long$Seed_Rain == 0),]


# 0 = no seeds/seedlings and 1 = seeds/seedlings
cover_rain_dis_long_unique   <- cover_rain_dis_long_unique   %>% 
 mutate(Aboveground = if_else(Aboveground > 0, 1, 0)) %>% 
 mutate(Seed_Rain = if_else(Seed_Rain > 0, 1, 0))

# Make a new column that adds both the seed rain and aboveground presence/absence columns. Then remove rows that = 2 (shares both species)
cover_rain_dis_long_unique$shared <- cover_rain_dis_long_unique$Aboveground + cover_rain_dis_long_unique$Seed_Rain
  
cover_rain_dis_long_unique2 <- cover_rain_dis_long_unique[!(cover_rain_dis_long_unique$shared == 2),]

cover_rain_dis_long_unique2 <- cover_rain_dis_long_unique2[, -6]


# Count unique taxa per site and transect for the seed bank compared to the seed rain

Unique_species_aboveground_com_rain <- cover_rain_dis_long_unique2[, -5]


Unique_species_aboveground_com_rain <- Unique_species_aboveground_com_rain %>% 
  filter(Aboveground != 0 ) %>% 
  group_by(Site, Transect) %>% 
  summarize(New_Aboveground = length(unique(SPP6)))

## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.

# Count unique taxa per site and transect for the seed rain compared to the seed bank

Unique_species_rain_com_aboveground <- cover_rain_dis_long_unique2[, -4]


Unique_species_rain_com_aboveground <- Unique_species_rain_com_aboveground %>% 
  filter(Seed_Rain != 0 ) %>% 
  group_by(Site, Transect) %>% 
  summarize(New_Rain = length(unique(SPP6)))

## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.

# Count the number of shared taxa between communities

cover_rain_dis_long_shared <- cover_rain_dis_long_unique %>% 
  filter(shared == 2) %>% 
  group_by(Site, Transect) %>% 
  summarize(Shared = length(unique(SPP6)))

## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.

## Merge them together

Species_counts_cover_rain <- full_join(Unique_species_aboveground_com_rain, Unique_species_rain_com_aboveground) %>% 
    mutate_all(~replace(., is.na(.), 0)) 

## Joining with `by = join_by(Site, Transect)`
## `mutate_all()` ignored the following grouping variables:

## Warning: There were 4 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `Transect = (structure(function (..., .x = ..1, .y = ..2, . =
##   ..1) ...`.
## ℹ In group 1: `Site = "PFCA 1"`.
## Caused by warning in `[<-.factor`:
## ! invalid factor level, NA generated
## ℹ Run ]8;;ide:run:dplyr::last_dplyr_warnings()dplyr::last_dplyr_warnings()]8;; to see the 3 remaining warnings.

Species_counts_cover_rain  <- full_join(Species_counts_cover_rain, cover_rain_dis_long_shared) %>% 
    mutate_all(~replace(., is.na(.), 0)) 

## Joining with `by = join_by(Site, Transect)`
## `mutate_all()` ignored the following grouping variables:

## Warning: There were 4 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `Transect = (structure(function (..., .x = ..1, .y = ..2, . =
##   ..1) ...`.
## ℹ In group 1: `Site = "PFCA 1"`.
## Caused by warning in `[<-.factor`:
## ! invalid factor level, NA generated
## ℹ Run ]8;;ide:run:dplyr::last_dplyr_warnings()dplyr::last_dplyr_warnings()]8;; to see the 3 remaining warnings.

Species_counts_cover_rain$Tot_Richness <- Species_counts_cover_rain$New_Aboveground + Species_counts_cover_rain$New_Rain + Species_counts_cover_rain$Shared


# Count the number of species in the seed rain for each transect

Species_counts_cover_rain$Tot_Rain <- Species_counts_cover_rain$New_Rain + Species_counts_cover_rain$Shared


# Count the number of species in the seed bank for each transect

Species_counts_cover_rain$Tot_Aboveground <- Species_counts_cover_rain$New_Aboveground + Species_counts_cover_rain$Shared


## Immigrants

Species_counts_cover_rain$Prop_Immigrant <- Species_counts_cover_rain$New_Rain / Species_counts_cover_rain$Tot_Aboveground

favstats(Species_counts_cover_rain$Prop_Immigrant~Species_counts_cover_rain$Site)

##   Species_counts_cover_rain$Site       min        Q1    median        Q3
## 1                         PFCA 1 0.2325581 0.2454268 0.3009259 0.4000000
## 2                         PFCA 2 0.1621622 0.2231378 0.2889714 0.3224343
## 3                         PFCA 3 0.1081081 0.2181818 0.3042017 0.3495879
## 4                             TP 0.1851852 0.2162162 0.2592593 0.3103448
##         max      mean         sd  n missing
## 1 0.6190476 0.3420640 0.12881067 10       0
## 2 0.3555556 0.2738325 0.06337656 10       0
## 3 0.4736842 0.2832612 0.11509477 10       0
## 4 0.5161290 0.2849422 0.10421660  9       0

# No difference in proportion of immigrant species
## *Note* better modeled as a binomial regression?

Species_counts_cover_rain.mod <- lm(Prop_Immigrant~Site, data =Species_counts_cover_rain)
summary(Species_counts_cover_rain.mod)

## 
## Call:
## lm(formula = Prop_Immigrant ~ Site, data = Species_counts_cover_rain)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.17515 -0.07599 -0.01296  0.05072  0.27698 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.34206    0.03345  10.226  4.7e-12 ***
## SitePFCA 2  -0.06823    0.04730  -1.442    0.158    
## SitePFCA 3  -0.05880    0.04730  -1.243    0.222    
## SiteTP      -0.05712    0.04860  -1.175    0.248    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1058 on 35 degrees of freedom
## Multiple R-squared:  0.06862,    Adjusted R-squared:  -0.01121 
## F-statistic: 0.8596 on 3 and 35 DF,  p-value: 0.4711

Anova(Species_counts_cover_rain.mod, type = "III")

## Anova Table (Type III tests)
## 
## Response: Prop_Immigrant
##              Sum Sq Df  F value    Pr(>F)    
## (Intercept) 1.17008  1 104.5809 4.703e-12 ***
## Site        0.02885  3   0.8596    0.4711    
## Residuals   0.39159 35                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ggplot(Species_counts_cover_rain, aes(x = Site, y = Prop_Immigrant)) +
  geom_point()+
  labs(x = "", y = "Proportion of immigrant species in seed rain", title = "Seed rain - Aboveground")

Cover vs. Seed Rain - Species-level and traits

By Site

# Make a new dataset to log the abundance variables

cover_rain_log <- cover_rain


# Do a log transformation (add the 0.1 to avoid infinite values)

cover_rain_log$Cover <- log(cover_rain_log$Cover+0.1)
cover_rain_log$Number_Seeds <- log(cover_rain_log$Number_Seeds+0.1)


# Get rid of infinite values

cover_rain_log[mapply(is.infinite, cover_rain_log)] <- 0

#cover_rain$Site <- factor(cover_rain$Site, levels = c("PFCA 1", "PFCA 2", "PFCA 3", "TP"), labels = c("Young", "Middle", "Old", "Remnant"))

ggplot(cover_rain_log, aes(x = Cover, y = Number_Seeds))+
  geom_point()+
  facet_wrap(vars(Site), nrow = 1, strip.position = "top")+
  geom_smooth(method = "lm") +
  labs(x = "Log(veg. cover)", y = "Log (no. seeds)")+
  theme(axis.text = element_text( size = 12),
        axis.text.x = element_text( size = 12),
        axis.title = element_text( size = 16),
        legend.position = "none",
        strip.text = element_text(size = 16))

## `geom_smooth()` using formula = 'y ~ x'

Same Plot

ggplot(cover_rain_log, aes(x = Cover, y = Number_Seeds, color = Site))+
  geom_point()+
  geom_smooth(method = "lm") +
  labs(y = "Log (Total Captured Seed Rain )", x = "Log (Total Cover )", title = "Seed rain vs cover")+
  theme_classic() 

## `geom_smooth()` using formula = 'y ~ x'

Residual plots

# Remove everything that we didn't see both cover and seeds for at a transect
## Removing the "rare" occurrences (cover has to be greater than 1 and number of seeds greater than 10) removes the streaking apparent in the residuals


cover_rain_filtered <- cover_rain %>% 
  filter(Cover > 0) %>% 
  filter(Number_Seeds > 0)

Fitting the model

Note I tried other distributional families (gamma, inverse guassian, tweedie, negative binomial) and logging both variables and using a Normal distribution seems to be the best option when looking at the residual vs. fits plot. You do get weird streaks in the residuals though because there were so many instances of us seeing 1 or 2 seeds/cover.

# Fit a linear mixed-model
mod <- lmer(log10(Number_Seeds)~log10(Cover)+(1|Site), data = cover_rain_filtered )

plot(mod)

Histograms of the predictor and response variables after logging show why we get the streaks. They still aren’t particularly Normal after logging.

p <-ggplot(cover_rain_filtered , aes(x=(log10(Number_Seeds)))) + 
  geom_histogram(color="black", fill="white")+
  labs(x = "Log(no. seeds)", y = "Frequency")
p

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

q <-ggplot(cover_rain_filtered , aes(x=log10(Cover))) + 
  geom_histogram(color="black", fill="white")+
  labs(x = "Log(total cover)", y = "Frequency")
q 

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Plotting the species and their associated residuals shows which species were over and under represented in the seed rain based on local vegetative cover. It does not show instances where we caught immigrating seeds (no local cover) or species with unsuccessful reproduction (no captured seeds).

resid.cover.rain <- residuals(mod)
SPP6 <- cover_rain_filtered $SPP6
Site <- cover_rain_filtered $Site

resid.species.all.cover.rain  <- data.frame(SPP6, resid.cover.rain, Site)
resid.specieslall.cover.rain  <- resid.species.all.cover.rain  %>% 
  arrange(desc(resid.cover.rain))

#Of the species captured in both the seed rain and cover, make a plot showing which species were present in the seed rain in greater numbers than expected. 


resid.species.all.cover.rain  %>% 
  mutate(SPP6 = fct_reorder(SPP6, resid.cover.rain, .fun='median')) %>% 
  ggplot(aes(x = SPP6, y = resid.cover.rain, color = Site)) +
  labs(x = "Species", y = "Residual", title = "Seed rain vs. Cover")+
  geom_point() +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 70, hjust = 1))+
  geom_hline(yintercept = 0, linetype = "dashed", color = "red", size = 1)+
  ylim(-5,5)+
  theme(axis.text = element_text( size = 12),
        axis.text.x = element_text( size = 5),
        axis.title = element_text( size = 16),
        strip.text = element_text(size = 16)) +
  annotate("text", x = 22, y = 3, label = "Overrepresented in seed rain")+annotate("text", x = 75, y = -4, label = "Underrepresented in seed rain")

Here is a plot showing which species were present in the seed rain in greater numbers than expected based on local cover. Bars represent 95% confidence intervals. Absent bars for a species indicates it was only present in the cover and seed rain in one transect.

Remember that this figure only shows species that we were able to captured both in the cover and the seed rain at a transect. Does not show anything we caught only seeds or only cover at.

The figure does line up very well with the results of Maya’s seed predation study. Some of the species favored by rodents were Helianthus mollis, Sporobolus heterolepis, Desmodium canadense, Liatris aspera, and Baptisia bracteata. Based on the figure, these species or congeners also had far fewer seeds that you would expect based on local cover. Other species on this end of the spectrum are also well known to be favored by birds (e.g., American goldfinch and bobwhite). For example, Cirsium altissimum, Silphium integrifolium, Echinacea pallida, Strophostyles leiosperma, and Ambrosia artemisiifolia. Baptisia bracteata and Baptisia alba are also known to suffer heavy pre-dispersal predation from insects at Tucker Prairie (Haddock et al. 1982).

Lysimachia lanceolata is interesting since it has a tight relationship with a rare genera of bees, Macropis. It is unclear how dependent this species is on pollination by Macropis for reproduction.

Size, both plant and seed, are also likely important factors. Based on Moles et al 2004, small plants also tend to produce many small seeds (e.g., Juncus) while large plants (e.g., Baptisia) produce fewer, large seeds with tradeoffs existing between seed size and juvenile survival.

Primary dispersal mode is likely also important here - many species high in cover but low in captured seeds have fleshy fruits. For example, Potentilla simplex, Rhus glabra, Rhus copallinum, Rosa spp., and Fragraria virginica. I believe the only fleshy fruited species we captured was Rubus spp.

Life history strategy is also important, with annuals being grouped on the overrepresented end of the spectrum.

resid.species.all.mean.cover.rain  <- resid.species.all.cover.rain  %>% 
  group_by(SPP6) %>% 
  summarize(mean.resid = mean(resid.cover.rain ),
            sd.resid = sd(resid.cover.rain),
            n.resid = n()) %>% 
  mutate(se.resid = sd.resid / sqrt(n.resid),
            lower.ci.resid = mean.resid - 1.96*((sd.resid/sqrt(n.resid))),
            upper.ci.resid = mean.resid + 1.96*((sd.resid/sqrt(n.resid))))

resid.species.all.mean.cover.rain  %>% 
  mutate(SPP6 = fct_reorder(SPP6, mean.resid, .fun='median')) %>% 
  ggplot(aes(x = SPP6, y = mean.resid)) +
  xlab("Species")+
  ylab("Residual") +
  geom_point() +
  geom_errorbar(aes(ymin = lower.ci.resid, ymax = upper.ci.resid), width = 0.2, position = position_dodge(0.05))+
  theme_classic() +
  theme(axis.text.x = element_text(angle = 70, hjust = 1))+
  geom_hline(yintercept = 0, linetype = "dashed", color = "red", size = 1) + 
  ylim(-5,5) +
  theme(axis.text = element_text( size = 10),
        axis.text.x = element_text( size = 5),
        axis.title = element_text( size = 16),
        strip.text = element_text(size = 16))+
  annotate("text", x = 22, y = 1.75, label = "Overrepresented in local seed rain")+
  annotate("text", x = 47, y = -1, label = "Underrepresented in local seed rain")

Mass_coverage <- left_join(resid.species.all.mean.cover.rain, traits)

## Joining with `by = join_by(SPP6)`

Mass_coverage  <- Mass_coverage  %>% 
  mutate(Seed_Mass = if_else(Mean_1_Seed_Mass_g > 0, 1, 0)) %>% 
  mutate_all(~replace(., is.na(.), 0))

pane1 <- Mass_coverage   %>% 
  mutate(SPP6 = fct_reorder(SPP6, mean.resid, .fun='median')) %>% 
  ggplot(aes(x = SPP6, y = mean.resid, color = as.factor(Seed_Mass))) +
  xlab("Species")+
  ylab("Residual") +
  geom_point() +
  scale_color_manual(values = c("black", "blue"))+
  geom_errorbar(aes(ymin = lower.ci.resid, ymax = upper.ci.resid), width = 0.2, position = position_dodge(0.05))+
  theme_classic() +
  theme(axis.text.x = element_text(angle = 70, hjust = 1))+
  geom_hline(yintercept = 0, linetype = "dashed", color = "red", size = 1) + 
  ylim(-5,5) +
  theme(axis.text = element_text( size = 10),
        axis.text.x = element_text( size = 5),
        axis.title = element_text( size = 16),
        strip.text = element_text(size = 16))+
  annotate("text", x = 22, y = 3, label = "Overrepresented in local seed rain")+
  annotate("text", x = 75, y = -4, label = "Underrepresented in local seed rain")
  
pane1

—-

Seed Bank and Vegetative Cover Exploratory Analysis

cover_bank_sum <- cover_bank %>% 
  group_by(SPP6) %>% 
  summarize(tot.cover = sum(Cover), tot.seeds = sum(Number_Seedlings))

Germinated but didn’t observed in the cover

• Abutilon theophrasti
• Amaranthus tuberculatus
• Cardamine hirsuta
• Cardamine parviflora
• Digitaria sanguinalis
• Dipsacus laciniatus
• Euphorbia humistrata
• Hordeum pusillum
• Ipomoea hederacea
• Mollugo verticillata
• Panicum capillare
• Panicum dichotomiflorum
• Persicaria longiseta
• Platanus occidentalis (likely a GH contaminate)
• Plantago pusilla
• Poa chapmaniana
• Populus deltoides (likely a GH contaminate)
• Portulaca oleracea
• Setaria glauca
• Toxicodendron radicans
• Verbena hastata
• Veronica peregrina
• Verbascum thapsus

Almost everything we were able to germinate that was not in the local cover were weedy annual species. Some things like the Cardamine might have been in the cover but have such an early phenology that we missed them. The two tree species might be greenhouse contaminates since the roof was partially open and these trees were right next door.

Observed in the cover but didn’t germinate

Too many to list (112 taxa)

Clearly, only a small fraction of local species survive the transition from cover -> seed rain -> to germinable seed bank.

Few perennial species were present in the germinable seed bank.

Dissimilarity analysis

## # A tibble: 2,025 × 5
## # Groups:   Site, Transect, SPP6 [1,725]
##    Site   Transect SPP6   Type        Abundance
##    <chr>  <fct>    <chr>  <chr>           <dbl>
##  1 PFCA 1 1        ACAVIR Aboveground      11.5
##  2 PFCA 1 1        ALOCAR Aboveground       1  
##  3 PFCA 1 1        AMBART Aboveground       6  
##  4 PFCA 1 1        BARVUL Aboveground       3.5
##  5 PFCA 1 1        BROJAP Aboveground       1  
##  6 PFCA 1 1        CERSPP Aboveground       1  
##  7 PFCA 1 1        CHAFAS Aboveground     388  
##  8 PFCA 1 1        CHAFAS Seed_Bank        11  
##  9 PFCA 1 1        CONCAN Aboveground       1  
## 10 PFCA 1 1        DESILL Aboveground       1  
## # … with 2,015 more rows

# Code obtained and modified from Lauren and Anu's seed bank vs. vegetation paper
## Only included Bray-Curtis and Jaccard dissimilarity matrices since the others produced similar results to these two

##### run dissimilarity loops
cover_bank_dis_results <-matrix(nrow=0,ncol=4)
sites <-as.vector(unique(cover_bank_dis$Site))


for(s in 1:length(sites)){
  temp <-subset(cover_bank_dis, Site==sites[s])
  temp$SPP6 <- factor(as.character(temp$SPP6))
  
  transects <-as.vector(unique(temp$Transect))
  
  #for(b in 1:length(blocks)){
    #temp_b<-subset(temp, block==blocks[b])

    #plots<-as.vector(unique(temp_b$plot))
   
      for(t in 1:length(transects)){ 
        
        temp_t <- subset(temp, Transect==transects[t])

        #wide form - with a row for above and belowground
        plot_cast<-dcast(temp_t[,-c(1:2)],Type ~ SPP6, value.var="Abundance", fun.aggregate=mean,fill=0)
        
        #relative abundances
        trans_tot <- cbind(plot_cast$Type,decostand(plot_cast[,-1],"tot"))
        names(trans_tot)[1]<-"trt"
        
        trans_pa <- cbind(plot_cast$Type,decostand(plot_cast[,-1],"pa"))
        names(trans_pa)[1]<-"trt"
        
        
        distance_bray<-vegdist(trans_tot[,-1], method = "bray")
        distance_jaccard <- vegdist(trans_pa[,-1], method = "jaccard")

        

        new.row<-c(sites[s], transects[t],distance_bray[1],distance_jaccard[1])
        cover_bank_dis_results  <-rbind( cover_bank_dis_results, new.row)
    }
  print(s/length(sites))
}

## [1] 0.25
## [1] 0.5
## [1] 0.75
## [1] 1

row.names(cover_bank_dis_results )<-NULL
cover_bank_dis_results <- data.frame(cover_bank_dis_results )

names(cover_bank_dis_results) <- c("Site", "Transect", "dissimilarity_bray", "dissimilarity_jaccard")

## 
## Call:
## lm(formula = dissimilarity_bray ~ Site, data = cover_bank_dis_results)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.220753 -0.058252  0.008957  0.061361  0.165558 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.55757    0.02931  19.023  < 2e-16 ***
## SitePFCA 2   0.26498    0.04145   6.393 2.08e-07 ***
## SitePFCA 3   0.31432    0.04145   7.583 5.73e-09 ***
## SiteTP       0.27515    0.04145   6.638 9.84e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09269 on 36 degrees of freedom
## Multiple R-squared:  0.6679, Adjusted R-squared:  0.6402 
## F-statistic: 24.13 on 3 and 36 DF,  p-value: 9.766e-09

## Anova Table (Type III tests)
## 
## Response: dissimilarity_bray
##              Sum Sq Df F value    Pr(>F)    
## (Intercept) 3.10879  1 361.886 < 2.2e-16 ***
## Site        0.62197  3  24.134 9.766e-09 ***
## Residuals   0.30926 36                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## Call:
## lm(formula = dissimilarity_jaccard ~ Site, data = cover_bank_dis_results)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.108471 -0.050240 -0.004958  0.054810  0.133643 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.78253    0.02066  37.884   <2e-16 ***
## SitePFCA 2   0.04594    0.02921   1.573   0.1246    
## SitePFCA 3   0.06965    0.02921   2.384   0.0225 *  
## SiteTP       0.05257    0.02921   1.800   0.0803 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06532 on 36 degrees of freedom
## Multiple R-squared:  0.1474, Adjusted R-squared:  0.07636 
## F-statistic: 2.075 on 3 and 36 DF,  p-value: 0.1208

## Anova Table (Type III tests)
## 
## Response: dissimilarity_jaccard
##             Sum Sq Df   F value Pr(>F)    
## (Intercept) 6.1236  1 1435.2106 <2e-16 ***
## Site        0.0266  3    2.0747 0.1208    
## Residuals   0.1536 36                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Fig. - Bray-Curtis dissimilarity

## 95% confidence intervals 

emmeans::emmeans(Cover_bank.mod.bray, "Site", type = "response")

##  Site   emmean     SE df lower.CL upper.CL
##  PFCA 1  0.558 0.0293 36    0.498    0.617
##  PFCA 2  0.823 0.0293 36    0.763    0.882
##  PFCA 3  0.872 0.0293 36    0.812    0.931
##  TP      0.833 0.0293 36    0.773    0.892
## 
## Confidence level used: 0.95

## Make a new dataframe to hold the model estimates and 95% confidence interval bounds

response <- c(0.558, 0.823, 0.872, 0.833)
lower <- c(0.498, 0.763,  0.812, 0.773) 
upper <- c(0.617, 0.882, 0.931,  0.892)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")


Cover_bank_dissimilarity <- data.frame(site, response, lower, upper)

Fig. - Jaccard dissimilarity figure

## 95% confidence intervals 

emmeans::emmeans(Cover_bank.mod.jac, "Site", type = "response")

##  Site   emmean     SE df lower.CL upper.CL
##  PFCA 1  0.783 0.0207 36    0.741    0.824
##  PFCA 2  0.828 0.0207 36    0.787    0.870
##  PFCA 3  0.852 0.0207 36    0.810    0.894
##  TP      0.835 0.0207 36    0.793    0.877
## 
## Confidence level used: 0.95

## Make a new dataframe to hold the model estimates and 95% confidence interval bounds

response <- c(0.783, 0.828, 0.852, 0.835)
lower <- c(0.741, .787,  0.810, 0.793) 
upper <- c( 0.824, 0.870, 0.894,  0.877)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")


Cover_bank_jac_dissimilarity <- data.frame(site, response, lower, upper)

Cover vs. Seed Bank Figures

By Site

# Drop Control for now

cover_bank <- cover_bank[!(cover_bank$Site%in%"CONTROL"),]


# Make a new dataset to log the abundance variables

cover_bank_log <- cover_bank


# Do a log transformation (add the 0.1 to avoid infinite values)

cover_bank_log$Cover <- log(cover_bank_log$Cover+0.1)
cover_bank_log$Number_Seedlings <- log(cover_bank_log$Number_Seedlings+0.1)

# Get rid of infinite values

cover_bank_log[mapply(is.infinite, cover_bank_log)] <- 0

cover_bank_log$Site <- factor(cover_bank_log$Site, levels = c("PFCA 1", "PFCA 2", "PFCA 3", "TP", "CONTROL"), labels = c("Young", "Middle", "Old", "Remnant", "Control"))

ggplot(cover_bank_log, aes(x = Cover, y = Number_Seedlings))+
  geom_point()+
  facet_wrap(vars(Site), nrow = 1, strip.position = "top")+
  geom_smooth(method = "lm") +
  labs(x = "Log(veg. cover)", y = "Log (no. seedlings)", title = "Seed bank vs. Cover")+
  theme(axis.text = element_text( size = 14),
        axis.text.x = element_text( size = 16),
        axis.title = element_text( size = 18),
        legend.position = "none",
        strip.text = element_text(size = 16))

## `geom_smooth()` using formula = 'y ~ x'

Same Plot

ggplot(cover_bank_log, aes(x = Cover, y = Number_Seedlings, color = Site))+
  geom_point()+
  geom_smooth(method = "lm") +
  labs(y = "Log (Total Germinated Seeds )", x = "Log (Total Cover )", title = "Seed bank vs cover")+
  theme_classic() 

## `geom_smooth()` using formula = 'y ~ x'

Residual plots

Fitting the Model

# Remove everything that we didn't see both cover and seed bank for at a transect

cover_bank_filtered <- cover_bank %>% 
  filter(Cover > 0) %>% 
  filter(Number_Seedlings > 0)



#
mod1 <- lmer(log(Number_Seedlings)~log(Cover)+(1|Site), data = cover_bank_filtered )


plot(mod1)

Again we get the streaks in the residual plots.

p1<-ggplot(cover_bank_filtered , aes(x=log(Number_Seedlings))) + 
  geom_histogram(color="black", fill="white")
p1

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

q1 <-ggplot(cover_bank_filtered , aes(x=log(Cover))) + 
  geom_histogram(color="black", fill="white")
q1 

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

resid.cover.bank <- residuals(mod1)
#resid <- filtered.mod$residuals


SPP6 <- cover_bank_filtered $SPP6
Site <- cover_bank_filtered $Site

resid.species.all.cover.bank <- data.frame(SPP6, resid.cover.bank, Site)
resid.specieslall.cover.bank <- resid.species.all.cover.bank %>% 
  arrange(desc(resid.cover.bank))

#Of the species captured in both the seed bank and cover, make a plot showing which species were present in the seed bank in greater numbers than expected. 

resid.species.all.cover.bank %>% 
  mutate(SPP6 = fct_reorder(SPP6, resid.cover.bank, .fun='median')) %>% 
  ggplot(aes(x = SPP6, y = resid.cover.bank, color = Site)) +
  labs(x = "Species", y = "Residual", title = "Seed Bank vs. Cover")+
  geom_point() +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 70, hjust = 1))+
  geom_hline(yintercept = 0, linetype = "dashed", color = "red", size = 1) +
  theme(axis.text = element_text( size = 10),
        axis.text.x = element_text( size = 8),
        axis.title = element_text( size = 16),
        strip.text = element_text(size = 16)) +
  ylim(-4,4)+
  annotate("text", x = 12, y = 3, label = "Overrepresented in seed bank")+annotate("text", x = 34, y = -4, label = "Underrepresented in seed bank")

Looking at this residual plot, disturbance tolerant annual species tend to be over represented in the seed bank compared to the local cover. These species are most likely enriched in seed banks overtime and are persistent. For example, Digitaria ischaemum, Setaria faberi, Plantago virginica, and Juncus spp. Some over representation might be because of superb dispersal ability (Eupatorium serotinum and Agrostis hyemalis).

Life history is clearly important, with annuals being more present than perennials.

The seed bank and cover had the least amount of overlap with shared species.

resid.species.all.mean.cover.bank <- resid.species.all.cover.bank %>% 
  group_by(SPP6) %>% 
  summarize(mean.resid = mean(resid.cover.bank),
            sd.resid = sd(resid.cover.bank),
            n.resid = n()) %>% 
  mutate(se.resid = sd.resid / sqrt(n.resid),
            lower.ci.resid = mean.resid - 1.96*((sd.resid/sqrt(n.resid))),
            upper.ci.resid = mean.resid + 1.96*((sd.resid/sqrt(n.resid))))

resid.species.all.mean.cover.bank %>% 
  mutate(SPP6 = fct_reorder(SPP6, mean.resid, .fun='median')) %>% 
  ggplot(aes(x = SPP6, y = mean.resid)) +
  labs(x = "Species", y = "Residual", title = "Seed Bank vs. Cover")+
  geom_point() +
  geom_errorbar(aes(ymin = lower.ci.resid, ymax = upper.ci.resid), width = 0.2, position = position_dodge(0.05))+
  theme_classic() +
  theme(axis.text.x = element_text(angle = 70, hjust = 1))+
  geom_hline(yintercept = 0, linetype = "dashed", color = "red", size = 1) + theme(axis.text = element_text( size = 12),
        axis.text.x = element_text( size = 8),
        axis.title = element_text( size = 16),
        strip.text = element_text(size = 16))+
  ylim(-4,4)+
  annotate("text", x = 10, y = 3, label = "Overrepresented in seed bank")+annotate("text", x = 34, y = -4, label = "Underrepresented in seed bank")

Mass_coverage_cover.bank <- left_join(resid.species.all.mean.cover.bank, traits)

## Joining with `by = join_by(SPP6)`

Mass_coverage_cover.bank  <- Mass_coverage_cover.bank  %>% 
  mutate(Seed_Mass = if_else(Mean_1_Seed_Mass_g > 0, 1, 0)) %>% 
  mutate_all(~replace(., is.na(.), 0))

pane2 <- Mass_coverage_cover.bank   %>% 
  mutate(SPP6 = fct_reorder(SPP6, mean.resid, .fun='median')) %>% 
  ggplot(aes(x = SPP6, y = mean.resid, color = as.factor(Seed_Mass))) +
  xlab("Species")+
  ylab("Residual") +
  geom_point() +
  scale_color_manual(values = c("black", "blue"))+
  geom_errorbar(aes(ymin = lower.ci.resid, ymax = upper.ci.resid), width = 0.2, position = position_dodge(0.05))+
  theme_classic() +
  theme(axis.text.x = element_text(angle = 70, hjust = 1))+
  geom_hline(yintercept = 0, linetype = "dashed", color = "red", size = 1) + 
  ylim(-5,5) +
  theme(axis.text = element_text( size = 10),
        axis.text.x = element_text( size = 5),
        axis.title = element_text( size = 16),
        strip.text = element_text(size = 16))+
  annotate("text", x = 16, y = 3, label = "Overrepresented in local seed bank")+
  annotate("text", x = 32, y = -4, label = "Underrepresented in local seed bank")
  
pane2

—-

Seed Rain and Seed Bank Exploratory Analysis

rain_bank_sum <- rain_bank %>% 
  group_by(SPP6) %>% 
  summarize(tot.seeds = sum(Number_Seeds), tot.seedlings = sum(Number_Seedlings))

Germinated but didn’t catch in the seed rain

• Abutilon theophrasti
• Digitalis sanguinalis
• Dipsacus laciniatus
• Euphorbia humistrata
• Lepidium densiflorum
• Panicum dichotomiflorum
• Plantago pusilla
• Poa chapmaniana

Mostly ruderal annual species.

Caught in the seed rain but didn’t germinate

Too many to list (61 taxa)

Dissimilarity analysis

## # A tibble: 1,859 × 5
## # Groups:   Site, Transect, SPP6 [1,500]
##    Site   Transect SPP6   Type      Abundance
##    <chr>  <fct>    <chr>  <chr>         <dbl>
##  1 PFCA 1 1        ACAVIR Seed_Rain        22
##  2 PFCA 1 1        AMBART Seed_Rain         6
##  3 PFCA 1 1        ANAMIN Seed_Rain        11
##  4 PFCA 1 1        BOLAST Seed_Rain         1
##  5 PFCA 1 1        CERSPP Seed_Rain        14
##  6 PFCA 1 1        CHAFAS Seed_Bank        11
##  7 PFCA 1 1        CHAFAS Seed_Rain       276
##  8 PFCA 1 1        CHEALB Seed_Rain         1
##  9 PFCA 1 1        CONCAN Seed_Rain         5
## 10 PFCA 1 1        DIGISC Seed_Bank         1
## # … with 1,849 more rows

# Code obtained and modified from Lauren and Anu's seed bank vs. vegetation paper
## Only included Bray-Curtis and Jaccard dissimilarity matrices since the others produced similar results to these two

##### run dissimilarity loops
rain_bank_dis_results <-matrix(nrow=0,ncol=4)
sites <-as.vector(unique(rain_bank_dis$Site))


for(s in 1:length(sites)){
  temp <-subset(rain_bank_dis, Site==sites[s])
  temp$SPP6 <- factor(as.character(temp$SPP6))
  
  transects <-as.vector(unique(temp$Transect))
  
for(t in 1:length(transects)){ 
      
        temp_t <- subset(temp, Transect==transects[t])

        #wide form - with a row for above and belowground
        plot_cast<-dcast(temp_t[,-c(1:2)],Type ~ SPP6, value.var="Abundance", fun.aggregate=mean,fill=0)
        
        #relative abundances
        trans_tot <- cbind(plot_cast$Type,decostand(plot_cast[,-1],"tot"))
        names(trans_tot)[1]<-"trt"
        
        trans_pa <- cbind(plot_cast$Type,decostand(plot_cast[,-1],"pa"))
        names(trans_pa)[1]<-"trt"
        
        
        distance_bray<-vegdist(trans_tot[,-1], method = "bray")
        distance_jaccard <- vegdist(trans_pa[,-1], method = "jaccard")

        

        new.row<-c(sites[s], transects[t],distance_bray[1],distance_jaccard[1])
        rain_bank_dis_results  <-rbind(rain_bank_dis_results, new.row)
    }
  print(s/length(sites))
}

## [1] 0.25
## [1] 0.5
## [1] 0.75
## [1] 1

row.names(rain_bank_dis_results)<-NULL
rain_bank_dis_results <- data.frame(rain_bank_dis_results)

names(rain_bank_dis_results) <- c("Site", "Transect", "dissimilarity_bray", "dissimilarity_jaccard")

## 
## Call:
## lm(formula = dissimilarity_bray ~ Site, data = rain_bank_dis_results)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.30162 -0.06832 -0.00551  0.07355  0.26438 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.55122    0.04061  13.575 9.94e-16 ***
## SitePFCA 2   0.24085    0.05743   4.194  0.00017 ***
## SitePFCA 3   0.18002    0.05743   3.135  0.00342 ** 
## SiteTP       0.13699    0.05743   2.386  0.02244 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1284 on 36 degrees of freedom
## Multiple R-squared:  0.3458, Adjusted R-squared:  0.2913 
## F-statistic: 6.344 on 3 and 36 DF,  p-value: 0.001449

## Anova Table (Type III tests)
## 
## Response: dissimilarity_bray
##              Sum Sq Df  F value    Pr(>F)    
## (Intercept) 3.03840  1 184.2741 9.941e-16 ***
## Site        0.31379  3   6.3436  0.001449 ** 
## Residuals   0.59359 36                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 
## Call:
## lm(formula = dissimilarity_jaccard ~ Site, data = rain_bank_dis_results)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.135047 -0.042571  0.002154  0.030196  0.128868 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.72479    0.02097  34.558  < 2e-16 ***
## SitePFCA 2   0.02556    0.02966   0.862  0.39458    
## SitePFCA 3   0.08522    0.02966   2.873  0.00677 ** 
## SiteTP       0.04026    0.02966   1.357  0.18317    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06632 on 36 degrees of freedom
## Multiple R-squared:  0.1949, Adjusted R-squared:  0.1278 
## F-statistic: 2.905 on 3 and 36 DF,  p-value: 0.04791

## Anova Table (Type III tests)
## 
## Response: dissimilarity_jaccard
##             Sum Sq Df   F value  Pr(>F)    
## (Intercept) 5.2532  1 1194.2473 < 2e-16 ***
## Site        0.0383  3    2.9051 0.04791 *  
## Residuals   0.1584 36                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Fig. - Bray-Curtis dissimilarity

## 95% confidence intervals 

emmeans::emmeans(Rain_bank.mod.bray, "Site", type = "response")

##  Site   emmean     SE df lower.CL upper.CL
##  PFCA 1  0.551 0.0406 36    0.469    0.634
##  PFCA 2  0.792 0.0406 36    0.710    0.874
##  PFCA 3  0.731 0.0406 36    0.649    0.814
##  TP      0.688 0.0406 36    0.606    0.771
## 
## Confidence level used: 0.95

## Make a new dataframe to hold the model estimates and 95% confidence interval bounds

response <- c(0.551, 0.792, 0.731, 0.688)
lower <- c(0.469, 0.710,  0.649, 0.606) 
upper <- c(0.634, 0.874, 0.814,  0.771)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")


Rain_bank_dissimilarity <- data.frame(site, response, lower, upper)

Fig. - Jaccard dissimilarity figure

## 95% confidence intervals 

emmeans::emmeans(Rain_bank.mod.jac, "Site", type = "response")

##  Site   emmean    SE df lower.CL upper.CL
##  PFCA 1  0.725 0.021 36    0.682    0.767
##  PFCA 2  0.750 0.021 36    0.708    0.793
##  PFCA 3  0.810 0.021 36    0.767    0.853
##  TP      0.765 0.021 36    0.723    0.808
## 
## Confidence level used: 0.95

## Make a new dataframe to hold the model estimates and 95% confidence interval bounds

response <- c(0.725, 0.750, 0.810, 0.765)
lower <- c(0.682, 0.708,  0.767, 0.723) 
upper <- c(0.767, 0.793, 0.853,  0.808)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")


Rain_bank_jac_dissimilarity <- data.frame(site, response, lower, upper)

Immigrants

rain_bank_dis_wide <- rain_bank_dis %>% 
  spread(SPP6, Abundance) %>% 
  mutate_all(~replace(., is.na(.), 0)) 

## `mutate_all()` ignored the following grouping variables:
## • Columns `Site`, `Transect`
## ℹ Use `mutate_at(df, vars(-group_cols()), myoperation)` to silence the message.

rain_bank_dis_long <- rain_bank_dis_wide %>% 
  gather(key = "SPP6", value = "Abundance", 4:144) %>% 
  spread(key = Type, Abundance) 

# Get rid of columns that both have zeros
rain_bank_dis_long_unique <- rain_bank_dis_long[!(rain_bank_dis_long$Seed_Bank == 0 & rain_bank_dis_long$Seed_Rain == 0),]


# 0 = no seeds/seedlings and 1 = seeds/seedlings
rain_bank_dis_long_unique   <- rain_bank_dis_long_unique   %>% 
  mutate(Seed_Bank = if_else(Seed_Bank > 0, 1, 0)) %>% 
   mutate(Seed_Rain = if_else(Seed_Rain > 0, 1, 0))

# Make a new column that adds both the seed rain and seed bank presence/absence columns. Then remove rows that = 2 (shares both species)
rain_bank_dis_long_unique$shared <- rain_bank_dis_long_unique$Seed_Bank + rain_bank_dis_long_unique$Seed_Rain
  
rain_bank_dis_long_unique2 <- rain_bank_dis_long_unique[!(rain_bank_dis_long_unique$shared == 2),]

rain_bank_dis_long_unique2 <- rain_bank_dis_long_unique2[, -6]


# Count unique taxa per site and transect for the seed bank compared to the seed rain

Unique_species_bank_com_rain <- rain_bank_dis_long_unique2[, -5]


Unique_species_bank_com_rain <- Unique_species_bank_com_rain %>% 
  filter(Seed_Bank != 0 ) %>% group_by(Site, Transect) %>% 
  summarize(New_Bank = length(unique(SPP6)))

## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.

# Count unique taxa per site and transect for the seed rain compared to the seed bank

Unique_species_rain_com_bank <- rain_bank_dis_long_unique2[, -4]


Unique_species_rain_com_bank <- Unique_species_rain_com_bank %>% 
  filter(Seed_Rain != 0 ) %>% 
  group_by(Site, Transect) %>% 
  summarize(New_Rain = length(unique(SPP6)))

## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.

# Count the number of shared taxa between communities

rain_bank_dis_long_shared <- rain_bank_dis_long_unique %>% 
  filter(shared == 2) %>% 
  group_by(Site, Transect) %>% 
  summarize(Shared = length(unique(SPP6)))

## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.

## Merge them together

Species_counts_rain_bank <- full_join(Unique_species_bank_com_rain, Unique_species_rain_com_bank) %>% 
    mutate_all(~replace(., is.na(.), 0)) 

## Joining with `by = join_by(Site, Transect)`
## `mutate_all()` ignored the following grouping variables:

## Warning: There were 4 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `Transect = (structure(function (..., .x = ..1, .y = ..2, . =
##   ..1) ...`.
## ℹ In group 1: `Site = "PFCA 1"`.
## Caused by warning in `[<-.factor`:
## ! invalid factor level, NA generated
## ℹ Run ]8;;ide:run:dplyr::last_dplyr_warnings()dplyr::last_dplyr_warnings()]8;; to see the 3 remaining warnings.

Species_counts_rain_bank  <- full_join(Species_counts_rain_bank, rain_bank_dis_long_shared) %>% 
    mutate_all(~replace(., is.na(.), 0)) 

## Joining with `by = join_by(Site, Transect)`
## `mutate_all()` ignored the following grouping variables:

## Warning: There were 4 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `Transect = (structure(function (..., .x = ..1, .y = ..2, . =
##   ..1) ...`.
## ℹ In group 1: `Site = "PFCA 1"`.
## Caused by warning in `[<-.factor`:
## ! invalid factor level, NA generated
## ℹ Run ]8;;ide:run:dplyr::last_dplyr_warnings()dplyr::last_dplyr_warnings()]8;; to see the 3 remaining warnings.

Species_counts_rain_bank$Tot_Richness <- Species_counts_rain_bank$New_Bank + Species_counts_rain_bank$New_Rain + Species_counts_rain_bank$Shared


# Count the number of species in the seed rain for each transect

Species_counts_rain_bank$Tot_Rain <- Species_counts_rain_bank$New_Rain + Species_counts_rain_bank$Shared


# Count the number of species in the seed bank for each transect

Species_counts_rain_bank$Tot_Bank <- Species_counts_rain_bank$New_Bank + Species_counts_rain_bank$Shared


## Immigrants

Species_counts_rain_bank$Prop_Immigrant <- Species_counts_rain_bank$New_Bank / Species_counts_rain_bank$Tot_Rain

ggplot(Species_counts_rain_bank, aes(x = Site, y = Prop_Immigrant)) +
  geom_point()

Seed Rain vs. Seed Bank Figures

By Site

# Drop the control for now

rain_bank <- rain_bank[!(rain_bank$Site%in%"CONTROL"),]

# Make a new dataset to log the abundance variables

rain_bank_log <- rain_bank


# Do a log transformation (add the 0.1 to avoid infinite values)

rain_bank_log$Number_Seedlings <- log(rain_bank_log$Number_Seedlings+0.1)
rain_bank_log$Number_Seeds <- log(rain_bank_log$Number_Seeds+0.1)


# Get rid of infinite values

rain_bank_log[mapply(is.infinite, rain_bank_log)] <- 0

rain_bank_log$Site <- factor(rain_bank_log$Site, levels = c("PFCA 1", "PFCA 2", "PFCA 3", "TP", "CONTROL"), labels = c("Young", "Middle", "Old", "Remnant", "Control"))

ggplot(rain_bank_log, aes(x = Number_Seeds, y = Number_Seedlings))+
  geom_point()+
  facet_wrap(vars(Site), nrow = 1, strip.position = "top")+
  geom_smooth(method = "lm") +
  labs(x = "Log(no. seeds)", y = "Log (no. seedlings)")+
  theme(axis.text = element_text( size = 14),
        axis.text.x = element_text( size = 20),
        axis.title = element_text( size = 18),
        legend.position = "none",
        strip.text = element_text(size = 16))

## `geom_smooth()` using formula = 'y ~ x'

Same Plot

ggplot(rain_bank_log, aes(x = Number_Seeds, y = Number_Seedlings, color = Site))+
  geom_point()+
  geom_smooth(method = "lm") +
  labs(y = "Log (Total Captured Seed Rain )", x = "Log (Total Cover )", title = "Seed rain vs Seed Bank")+
  theme_classic() 

## `geom_smooth()` using formula = 'y ~ x'

Residual plots

# Remove everything that we didn't see both the seed rain and bank at a transect

rain_bank_filtered <- rain_bank %>% 
  filter(Number_Seedlings> 0) %>% 
  filter(Number_Seeds > 0)

#
mod3 <- lmer(log(Number_Seedlings)~log(Number_Seeds)+(1|Site), data = rain_bank_filtered)


plot(mod3)

Again we get the streaks in the residual vs. fitted values plot.

p2<-ggplot(rain_bank_filtered , aes(x=log(Number_Seedlings))) + 
  geom_histogram(color="black", fill="white")
p2

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

q2 <-ggplot(rain_bank_filtered , aes(x=log(Number_Seeds))) + 
  geom_histogram(color="black", fill="white")
q2

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

resid.rain.bank<- residuals(mod3)

SPP6 <- rain_bank_filtered$SPP6
Site <- rain_bank_filtered$Site

resid.species.all.rain.bank <- data.frame(SPP6, resid.rain.bank, Site)
resid.specieslall.rain.bank <- resid.species.all.rain.bank  %>% 
  arrange(desc(resid.rain.bank))

#Of the species captured in both the seed rain and bank, make a plot showing which species were present in the seed bank in greater numbers than expected. 

resid.species.all.rain.bank %>% 
  mutate(SPP6 = fct_reorder(SPP6, resid.rain.bank, .fun='median')) %>% 
  ggplot(aes(x = SPP6, y = resid.rain.bank, color = Site)) +
  labs(x = "Species", y = "Residual", title = "Seed Rain vs. Seed Bank") +
  geom_point() +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 70, hjust = 1))+
  geom_hline(yintercept = 0, linetype = "dashed", color = "red", size = 1)+
  theme(axis.text = element_text( size = 12),
        axis.text.x = element_text( size = 8),
        axis.title = element_text( size = 16),
        strip.text = element_text(size = 16))+
  ylim(-4,4)+
  annotate("text", x = 12, y = 3, label = "Overrepresented in seed bank")+annotate("text", x = 40, y = -4, label = "Underrepresented in seed bank")

resid.species.all.mean.rain.bank <- resid.species.all.rain.bank %>% 
  group_by(SPP6) %>% 
  summarize(mean.resid = mean(resid.rain.bank),
            sd.resid = sd(resid.rain.bank),
            n.resid = n()) %>% 
  mutate(se.resid = sd.resid / sqrt(n.resid),
            lower.ci.resid = mean.resid - 1.96*((sd.resid/sqrt(n.resid))),
            upper.ci.resid = mean.resid + 1.96*((sd.resid/sqrt(n.resid))))

Interesting that some species underrepresented in the seed rain based on local cover are present in germinable seed bank (LYSLAN, CIRALT, STRLEI, EUTGYM, RATPIN). Does this suggest high mortality when transitionig into the seed bank and low mortality once in the seed bank?

Again, most species overrepresented in the seed bank based on local seed input are weedy and annual species. Most perennials fall on the underrepresented side of the spectrum.

Many of the native ones also seem to be species that do well in restorations (Lespedeza capitata, Ratibida pinnata, Andropogon gerardii, Sorghastrum nutans, Coreopsis lanceolata, Tradescantia ohiensis, Penstemon digitalis, Schizachyrium scoparium, Solidago altissima, Symphytrichum spp., etc.). It’s probably no coincidence that these species that recruit from seed rain also do well in restorations that are broadcast seeded. I wonder if restoration practices are accidentally selecting for species that recruit well from seed (limited dispersal opportunities and reinforcement from what initially recruits).

I think the most surprising species that germinated is Lysimachia lanceolata, which is not present in the PFCA restorations.

resid.species.all.mean.rain.bank %>% 
  mutate(SPP6 = fct_reorder(SPP6, mean.resid, .fun='median')) %>% 
  ggplot(aes(x = SPP6, y = mean.resid))  +
  labs(x = "Species", y = "Residual", title = "Seed Rain vs. Seed Bank") +
  geom_point() +
  geom_errorbar(aes(ymin = lower.ci.resid, ymax = upper.ci.resid), width = 0.2, position = position_dodge(0.05))+
  theme_classic() +
  theme(axis.text.x = element_text(angle = 70, hjust = 1))+
  geom_hline(yintercept = 0, linetype = "dashed", color = "red", size = 1) + theme(axis.text = element_text( size = 12),
        axis.text.x = element_text( size = 8),
        axis.title = element_text( size = 16),
        strip.text = element_text(size = 16))+
  ylim(-4,4)+
  annotate("text", x = 12, y = 3, label = "Overrepresented in seed bank")+annotate("text", x = 40, y = -4, label = "Underrepresented in seed bank")

resid.species.all.mean.rain.bank %>% 
  mutate(SPP6 = fct_reorder(SPP6, mean.resid, .fun='median')) %>% 
  ggplot(aes(x = SPP6, y = mean.resid)) +
  labs(x = "Species", y = "Residual", title = "Seed Bank vs. Cover")+
  geom_point() +
  geom_errorbar(aes(ymin = lower.ci.resid, ymax = upper.ci.resid), width = 0.2, position = position_dodge(0.05))+
  theme_classic() +
  theme(axis.text.x = element_text(angle = 70, hjust = 1))+
  geom_hline(yintercept = 0, linetype = "dashed", color = "red", size = 1) + theme(axis.text = element_text( size = 12),
        axis.text.x = element_text( size = 8),
        axis.title = element_text( size = 16),
        strip.text = element_text(size = 16))+
  ylim(-4,4)+
  annotate("text", x = 10, y = 3, label = "Overrepresented in seed bank")+annotate("text", x = 34, y = -4, label = "Underrepresented in seed bank")

Mass_coverage_rain.bank <- left_join(resid.species.all.mean.rain.bank, traits)

## Joining with `by = join_by(SPP6)`

Mass_coverage_rain.bank  <- Mass_coverage_cover.bank  %>% 
  mutate(Seed_Mass = if_else(Mean_1_Seed_Mass_g > 0, 1, 0)) %>% 
  mutate_all(~replace(., is.na(.), 0))

pane3 <- Mass_coverage_rain.bank  %>% 
  mutate(SPP6 = fct_reorder(SPP6, mean.resid, .fun='median')) %>% 
  ggplot(aes(x = SPP6, y = mean.resid, color = as.factor(Seed_Mass))) +
  xlab("Species")+
  ylab("Residual") +
  geom_point() +
  scale_color_manual(values = c("black", "blue"))+
  geom_errorbar(aes(ymin = lower.ci.resid, ymax = upper.ci.resid), width = 0.2, position = position_dodge(0.05))+
  theme_classic() +
  theme(axis.text.x = element_text(angle = 70, hjust = 1))+
  geom_hline(yintercept = 0, linetype = "dashed", color = "red", size = 1) + 
  ylim(-5,5) +
  theme(axis.text = element_text( size = 10),
        axis.text.x = element_text( size = 5),
        axis.title = element_text( size = 16),
        strip.text = element_text(size = 16))+
  annotate("text", x = 16, y = 3, label = "Overrepresented in local seed rain")+
  annotate("text", x = 32, y = -4, label = "Underrepresented in local seed rain")
  
pane3

—- Multipaned figures